Exploitation of Pilot Signals for Blind Resilient Detection and Geo-Observable Estimation of Navigation Signals

ABSTRACT

A method and apparatus detects and estimates geo-observables of navigation signals employing civil formats with repeating baseband signal components, i.e., “pilot signals,” including true GNSS signals generated by satellite vehicles (SV’s) or ground beacons (pseudolites), and malicious GNSS signals, e.g., spoofers and repeaters. Multi-subband symbol-rate synchronous channelization can exploit the full substantive bandwidth of the GNSS signals with managed complexity in each subband. Spatial/polarization receivers can be provided to remove interference and geolocate non-GNSS jamming sources, as well as targeted GNSS spoofers that emulate GNSS signals. This can provide time-to-first-fix (TTFF) over much smaller time intervals than existing GNSS methods; can operate in the presence of signals with much wider disparity in received power than existing techniques; and can operate in the presence of arbitrary multipath

CROSS REFERENCE TO RELATED APPLICATIONS

This application is a Continuation of U.S. Pat. Appl. No. 17/479,658, filed Sep. 20, 2021, now U.S. Pat. No. 11,650,328; which is a Continuation of U.S. Pat. Appl. No. 16/701,144, filed on Dec. 2, 2019, now U.S. Pat. No. 11,125,888; which claims priority to Provisional Appl. No. 62/773,589, filed Nov. 30, 2018, and Provisional Appl. No. 62/773,605, filed Nov. 30, 2018; and U.S. Pat. Appl. No. 16/701,144 is a Continuation-in-Part of U.S. Pat. Appl. No. 15/731,417, filed Jun. 5, 2017, now U.S. Pat. No. 10,775,510; which claims priority to Provisional Appl. No. 62/392,623, filed Jun. 6, 2016, and Provisional Appl. No. 62/429,029, filed Dec. 1, 2016; all of which are hereby incorporated by reference in their entireties and all of which this application claims priority under at least 35 U.S.C. 120 and/or any other applicable provision in Title 35 of the United States Code.

FIELD OF THE INVENTION

Aspects of the disclosure relate generally to location determination, and more particularly to detecting or determining geo-observables of navigation signals transmitted by satellite vehicles (SV’s) in Global Navigation Satellite Systems (GNSS), airborne non-SV beacons (e.g., pseudolites), and ground-based beacons.

BACKGROUND

The American Global Positioning System (GPS) L1 coarse acquisition signal (commonly referred to as the “L1 C/A signal,” or more recently as the “L1 legacy signal”), originally developed as a means for aiding acquisition of the military-grade P (later encrypted P(Y)) code, has gained near-ubiquitous use as a means for navigation and carrier/timing synchronization in civil positioning, navigation, and timing (PNT) systems. This is due in part to the GPS network’s worldwide visibility; the publicly disclosed and well-known structure of the L1 C/A signal-in-space (SiS); and the ready availability of mature, low-cost hardware for receiving the L1 C/A signal.

The success of the GPS L1 C/A signal has resulted in incorporation of “civil” modes in all of the Global Navigation Satellite System (GNSS) signals developed to date, e.g., the GLONASS, Galileo, Beidou, and NAVIC/IRNSS signals. In addition, the need for additional precision to support civil PNT operation in urban canyons and autonomous navigation systems has spurred the development of the wideband civil GPS L5, Galileo E5B, and E6B/C modes, and the Japanese Quasi-Zenith Satellite System (QZSS), which is transmits a version of the civil GPS L5 signal from geosynchronous orbit. Moreover, similar civil signal modes have been proposed as GNSS-like ground and air beacons, such as the Locata’s LocataLite signal, which emulates the GPS C/A signal with a ×10 spreading rate in the 2.4-2.485 GHz ISM band. All of these civil modes possess a common modulation format, referred to herein as modulation-on-symbol direct-sequence spread spectrum (MOS-DSSS), in which a baseband symbol stream, typically (but not always) operating at a 1,000 sample-per-second (1 ksps) symbol rate, is modulated by a higher-rate ranging code that repeats every symbol.

Since its inception, the only method used to acquire, demodulate, and determine geo-observables of these navigation signals in PNT applications has been matched-filter despreading methods. These methods operate by correlating the ranging code (or a set of ranging codes that the transmitter may be using) against the received signal and adjusting the carrier and timing offset of that code (referred to herein as geo-observables of the received navigation signal) until it aligns with the code received from the GNSS transmitter. However, matched filter despreading possesses inherent limitations that can degrade its performance in many reception environments. Known limitations include the following:

-   Slow cold-start time-to-first-fix (TTFF), even in the most benign     environments. In a cold start scenario, the receiver does not know     the specific ranging codes used by any of the GNSS satellite     vehicles, nor does it know the timing phase or frequency offset of     those codes, hence it must search over code, timing, and frequency     to synchronize with a single satellite. Cold-start TTFF can be on     the order of minutes, and in addition can be quite power     consumptive. -   High sensitivity to multipath interference, including inadvertent     multipath induced by the specular reflectors in the vicinity of the     receiver, which can distort the ranging codes and degrade or destroy     correlation peaks, and synthetic multipath induced by malicious     repeaters. Inadvertent multipath is especially acute in urban and     indoor navigation and localization scenarios. -   High susceptibility to strong co-channel interference, including     near-far interference from local beacons, malicious non-targeted     spoofing signals, and malicious jamming signals. The GNSS signals     deployed to date only provide a despreading gain of 30-40 dB, based     on the chip-rate of the GNSS ranging codes, most of which is needed     to raise the GNSS signal above the background noise; hence the     signal is easily jammed by many types of interference. An additional     10-13 dB gain can be obtained by exploiting additional structure of     the baseband symbol stream. However, this further increases TTFF in     cold-start acquisition scenarios. -   High susceptibility to targeted or “aligned” covert spoofers that     can use knowledge of the ranging code from the satellite and rough     location of the receiver platform to duplicate and attempt to     supplant signals from those satellites. -   Reliance on ancillary anti-spoofing (AS) codes, e.g., the GPS L1     P(Y) code, to overcome jamming and spoofing events. This reliance     increases the cost of the receiver, limits AS capability to a subset     of military receivers (again, with inherently higher cost), and     introduces a host of physical security measures and restrictions to     prevent compromise of AS-capable receivers. -   Need to maintain a library of all of the GNSS ranging codes, or     seeds used to initiate the codes, in order to implement the     despreader. This can complicate the receiver architecture,     especially for “multi-GNSS receivers” that are designed to exploit     civil signals generated by multiple families of satellites.

Moreover, there is a perceived (but not real) need for matched-filter despreaders to receive and sample signals over wide bandwidths, constrained to a large integer multiple of the ranging code chip-rate, e.g., 1.023 million chips per second (Mcps) for the GPS C/A code, and 10.23 Mcps for the GPS L5 short code, in order to fully exploit the processing gain of the code, overcome aforementioned channel multipath effects, and (in military GPS receivers) exploit moderate and long codes overlaid on top of the short civil code. This perceived need is further complicated in multi-GNSS radios, which must accommodate wide ranges of chip rates and bandwidths.

SUMMARY

The following summarizes some aspects of the present disclosure to provide a basic understanding of the discussed technology. This summary is not an extensive overview of all contemplated features of the disclosure, and is intended neither to identify key or critical elements of all aspects of the disclosure nor to delineate the scope of any or all aspects of the disclosure. Its sole purpose is to present some concepts of one or more aspects of the disclosure in summary form as a prelude to the more detailed description that follows.

Disclosed aspects provide for detecting and obtaining geo-observables of navigation signals generated by satellite vehicles (SV’s) in Global Navigation Satellite Systems (GNSS), and generated by ground-based or airborne non-SV beacons, e.g., pseudolites. This includes effecting detection and geo-observable estimation in environments populated with many sources, including malicious sources intending to disrupt or subvert information provided by legitimate sources, for example, non-GNSS jammers, spoofers that emulate GNSS signals, and repeaters that can record and replay environments containing true GNSS signals. Furthermore, this includes effecting said detection and geo-observable estimation in environments subject to potentially severe differences in power levels between such signals, for example, “near-far interference” between remote GNSS signals and nearby ground or air based beacons, jammers, and spoofers.

Aspects of the disclosure can achieve these and other goals by exploiting the massive spectral redundancy of navigation signals employing modulation-on-symbol direct-sequence spread spectrum (MOS-DSSS) modulation formats, in which a baseband symbol sequence is spread by a ranging code that repeats every each baseband symbol. Novel and powerful methods for exploiting this signal structure are introduced. For example, aspects described herein can exploit the time, frequency and code diversity of those signals, to implement adaptive linear combining methods and linear-algebraic combiner adaptation algorithms to detect all of the MOS-DSSS navigation signals in the received environment at the maximum signal-to-interference-and-noise ratio (max-SINR) achievable by the combiners; estimate key geo-observables of those detected signals; identify detected malicious sources ahead of or during PNT acquisition and tracking operations; prevent those malicious sources from corrupting or subverting the PNT operation; and optionally alert the receiver to the presence of those sources, and/or geolocate those malicious sources. In aspects employing multifeed receivers wherein multiple spatial and/or polarization diverse antennas are coupled to multiple receivers that are receptive to the MOS-DSSS navigation signals, aspects of the disclosure can additionally exploit the additional spatial and/or polarization diversity of navigation and non-navigation (e.g., jamming) signals to further remove non-navigation signals from the received MOS-DSSS navigation signals, and to remove targeted spoofing signals that closely emulate the time, frequency, and code diversity of navigation signals the spoofers are intending to displace at a victim receiver.

Moreover, such aspects can perform all of these operations without knowledge of, a search over, or synchronization to the specific ranging codes employed by each navigation signal. As a consequence, the system described herein can enable faster time-to-first-fix (TTFF) than conventional navigation receivers that require a search over ranging code parameters, and additional robustness to multipath interference (including synthetic multipath induced by repeaters) that can impair that search. Disclosed aspects also provide enhanced compatibility with implementations employing software-defined radio (SDR) architectures, including general-purpose SDR (GPSDR) architectures that are not optimized to GNSS applications, further improving cost and flexibility of the system.

Aspects disclosed herein can adapt symbol-rate (e.g., not chip-rate) synchronous reception and channelization structure, using both purpose-built equipment designed specifically for GNSS bands and signals of interest, or general-purpose software-defined receiver (GPSDR) equipment, e.g., Ettus Research Universal Research Software Peripheral (USRP) radios such as the B200mini-series device, which can employ a much wider range of sampling rates and bandwidths, and which employ analog-to-digital converters (ADC) with 8 or more bits per in-phase (1) and quadrature (Q) rail, allowing implementation of sophisticated jammer mitigation techniques and operation in severe multipath environments. This is in contrast to GNSS-specific SDR (GNSS-SDR) equipment, e.g., the Maxim MAX2769 series of Universal GNSS Receivers, which are employ ADC’s with sampling rates preset to integer multiple of the GNSS chip-rate, e.g., 16.238 Msps (16×1.023 Msps), and with ADC precisions of 1-2 bits per I and Q rail, or 3 bits per I rail, thereby substantively worsening the receivers’ vulnerability to jamming measures and severe multipath.

Aspects disclosed herein can be configured to operate in devices, systems, and methods disclosed in nonprovisional patent application 15/731,417, “Blind Despreading of Civil GNSS Signals for Resilient PNT Applications” (the ‘417 NPA), incorporated herein by reference. Such aspects can be employed in a method, system, and/or apparatus for receiving, blindly despreading, and determining geo-observables, of true civil Global Navigation Satellite Systems (GNSS) navigation signals generated by any of the set of satellite vehicles and ground beacons, amongst false echoes and malicious GNSS signals from spoofers and repeaters; for identifying malicious GNSS signals, and preventing those signals from corrupting or capturing PNT tracking operations; and for geolocating malicious GNSS signals. Some aspects can be configured for spatial/polarization diverse receivers that remove non-GNSS jammers received by the system, as well as targeted GNSS spoofers that can otherwise emulate GNSS signals received at victim receivers.

Aspects of the disclosure are further described for exploiting navigation signal modulation formats in which a component of the baseband symbol sequence, referred to herein as a “pilot signal,” is periodically repeated. Novel and powerful methods for exploiting this signal structure is an additional focus of aspects disclosed herein. For example, means are described for exploiting the period or content of those pilot signals to effect “partially blind” signal detection and geo-observable estimation methods that are highly efficient relative to equivalent “fully blind” methods that exploit more general properties of the baseband symbol sequence; and that can perform this detection and geo-observable estimation at much lower received incident powers (RIP’s) than equivalent fully-blind methods.

Aspects of the disclosure are further described for exploiting navigation or symbol streams provided by third parties. It is known that fully-blind despreading methods introduce complexity inherent to the goals and needs of the method; specifically, the demodulation (symbol estimation) step. This step inherently requires a despreader output signal-to-noise ratio (SNR) that is high enough to determine the baseband symbols after integration over each symbol, or to demodulate the navigation sequence after integration over each navigation symbol. This introduces a fundamental “acquisition threshold” for any fully blind method, i.e., the received incident power (RIP) or despread signal-to-inference-and-noise ratio (SINR) at which navigation signals can be reliably demodulated by the methods, or for that matter for any conventional matched-filter despreading method. This also places a limit on the accuracy of any geo-observable estimation algorithm that relies on that demodulated navigation data, which will be degraded by errors in that data.

However, in many use scenarios, baseband symbol or navigation sequences can be obtained from third-party sources. For example, in some use scenarios, PNT-capable platforms may possess GNSS reception devices that can receive, despread, and demodulate navigation data on their own in non-challenged environments. For example, recent versions of Android provide raw GPS navigation bit-streams along with GNSS measurement results. In other use scenarios, a PNT-capable or noncapable platform may possess communication links to other PNT-capable platforms or network infrastructure capable of receiving, despreading, and demodulating navigation signals. For example CORS (Continuously-Operating Reference Stations) decode navigation data and make it available in a decoded format that can be obtained over the Internet.

This third-party data can be used to implement “partially-blind” methods that can achieve the blind signal separation, jammer/spoofer resilience, near-far interference resilience, and multipath robustness without the need for, or sensitivity to, demodulating the received navigation data. Among other advantages, these approaches can develop resilient PNT (RPNT) analytics over integration times that are substantively (in theory, arbitrarily) longer than a single baseband or navigation symbol, allowing accurate signal geo-observable estimation well below the acquisition threshold required for both any fully-blind despreading method. In many use scenarios, this performance can furthermore be provided at substantively lower complexity than fully-blind methods that require demodulation of the navigation sequence.

Aspects disclosed herein also allow the use of multi-subband architectures in which symbol-rate synchronous channelizers are employed over multiple subbands within the GNSS signal band, or symbol-rate synchronous channelizer bins covering the GNSS signal bandwidth are organized into subbands within the GNSS signal band. In these aspects, linear combiners, and interference-excising linear-algebraic combiner adaptation algorithms, are employed to develop detection features and signal estimates within each subband metrics, which are then combined across subbands. The resultant method can exploit the full bandwidth of the GNSS signal, with a system complexity that scales linearly with the bandwidth exploited (number of subbands), and can constrain the dimensionality (and therefore complexity) of linear combiner adaptation operations employed in each subband to just that needed to effectively excise interference encountered within that subband. The multi-subband approach also provides a means for naturally excising non-tonal narrowband jamming signals, by selectively and adaptively suppressing subbands containing strong non-GNSS interference.

In one aspect, a low-cost single-feed reception system and multi-subband symbol-rate synchronous channelizer is presented that can develop resilient positioning, navigation, and timing (RPNT) analytics for civil modes of Global Navigation Satellite System (GNSS) signals generated by satellite vehicles (SV’s), civil modes of navigation signals generated by ground and air beacons generated by pseudolites, in the presence of large numbers of malicious tonal and narrowband jammers; malicious signals attempting to spoof those signals; and time and frequency selective multipath induced by natural reflections or malicious repeaters. The system can operate without prior knowledge of, or a search over, the ranging codes for those signals, and even if those signals are received at widely separated power levels; and perform this differentiation based only on the time, frequency, and code diversity between those signals.

In one aspect, the multi-subband channelization architecture is extended to spatial and/or polarization diverse multi-feed reception systems, in which each subband employs linear combiners, and linear-algebraic combiner adaptation algorithms, which combine data over spectral and spatial/polarization diversity dimensions to separate GNSS signals and excise jamming interference, and then combine detection features from each subband over the full bandwidth of the GNSS signal. Such systems can exploit spatial and/or polarization diversity between the true navigation signals and the malicious spoofers and repeaters to further blindly detect and differentiate the navigation signals from wideband jammers, and from signals generated by closely aligned repeaters and targeted spoofers that are attempting to match the time, frequency, and code diversity of the navigation signals. As an added advantage, the dimension of each subband combiner can be held constant as the number of spatial/polarization diversity dimensions grows, resulting in a linear growth in complexity as spatial/polarization diversity is increased.

Aspects of the disclosure described herein include a system having identifying means comprising at least one antenna, coupled to at least one receiver receptive to energy in a band containing at least one navigation signal with an MOS-DSSS modulation format. The identifying means can include a receiver, for example, dedicated or software-defined radio resources, and a digital signal processor (DSP) (e.g., DSP FPGA, GPU, or the like, with associated storage capability), or combinations of DSP hardware configured for implementing software or firmware installed on that DSP hardware. In these aspects, reception operations can comprise:

-   receiving and amplifying at least one signal-in-space (SiS)     containing at least a subband of the at least one navigation signal,     comprising at least one signal with a MOS-DSSS modulation format     comprising a baseband signal with symbol rate f_(sym), spread by a     ranging code that repeats every baseband symbol; -   frequency-shifting the at least one SiS to an intermediate or     complex baseband frequency, e.g., using a direct-conversion mixer     with LO frequency f_(R) and appropriate filtering operations; and -   sampling the frequency-shifted SiS using an analog-to-digital     convertor (ADC) with a sampling rate f_(ADC) that is equal to an     integer multiple of the baseband symbol rate f_(sym), i.e.,     f_(ADC)=M_(ADC)f_(sym′) where M_(ADC) is a (typically large)     integer.

In these aspects, the DSP can perform the following operations:

-   capturing at least one sampled and frequency-shifted SiS over a     block of received data samples, and recording the approximate     reception time and reception frequency of that data block, referred     to herein as a “time/frequency stamp”; -   channelizing the sampled received data block using at least one     symbol-rate synchronous channelizer to provide at least K_(chn)     frequency channels with output rate f_(sym′) each with frequency     extent lying within the active bandwidth of the MOS-DSSS navigation     signal; -   organizing the at least K_(chn) frequency channels into K_(sub)     subsets of frequency channels, referred to herein as subbands, each     comprising M_(sub) frequency channels; -   in multifeed receivers, wherein multiple spatial and/or polarization     diverse antennas are coupled to M_(feed) receivers receptive to     energy in a band containing at least one navigation signal with an     MOS-DSSS modulation format, channelizing each feed into K_(chn)     frequency channels with output rate f_(sym), organizing those     frequency channels into K_(sub) subbands each comprising     M_(sub)frequency channels, and combining those subbands across     receiver feeds to form a M_(DoF) _(×1) vector signal with symbol     rate f_(sym) for each subband, where M_(DoF) = M_(sub) M_(feed)is     referred to herein as the degrees of freedom (DoF’s) of the subband,     and where M_(DoF) = M_(sub) for the single-feed receiver (M_(feed) =     1); -   collecting the K_(sub) subbands of M_(DoF) ×1 data vectors over     N_(sym) symbols, and adding a time and frequency stamp those     collections of symbols, to form a time/frequency stamped     multi-subband output data block; -   obtaining third-party baseband symbol data or navigation data from     time channels and frequency bands covering the time/frequency stamp,     including uncertainty between the receiver LO, clock rate, and/or     clock timing phase and universal time coordinates (UTC), using     external resources also available to the system; and -   processing the time/frequency stamped multi-subband output data     block to obtain RPNT analytics, using blind adaptive processing     operations responsive to at least one navigation signal received by     the system, and using the third-party baseband symbol data or     navigation data.

In some aspects of the disclosure, the sampled receive data block, or the channelized received data block, along with its time/frequency stamp, might be transported to additional DSPs for processing to obtain the RPNT analytics. These additional DSPs may be onboard the platform containing the reception means, or in a different location entirely, including infrastructure equipment, i.e., “in the Cloud,” if the platform possesses communications means allowing those data blocks to be transported from the platform.

In some aspects, the reception operations are instantiated using purpose-built hardware and DSP hardware. In other aspects, the reception operations are instantiated using general-purpose software-defined radio (GPSDR) and DSP hardware. In either aspect, the receiver clock may be adjusted to remove error from UTC, once a PNT solution is obtained for the receiver. Alternately, the ADC output data may be frequency-sampled and resampled to conform to UTC using algorithms implemented in DSP hardware, once that solution is obtained.

In the various aspects, the ADC sampling rate need not have any relation to the chip-rate of the navigation signal; nor need it have any relation to the bandwidth of the navigation signal. In particular, it can have a value that is easily instantiated in low-cost, general-purpose, receiver hardware; and it can have a value that is substantively less than the bandwidth of the navigation signal.

In one aspect of the disclosure, the symbol-rate synchronous subband generation is performed by first applying an M_(ADC)-sample fast Fourier transform (FFT) to each S/P output data vector on each of M_(feed) receiver feeds, and then selecting K_(FFT) output bins lying within the active bandwidth of the MOS-DSSS navigation signal, to create K_(chn) frequency channels each with output rate f_(sym) and dimensionality M_(feed) ×1. These M_(feed) ×1 dimensional frequency channels are then organized into K_(sub) subbands with M_(sub) channels per subband, creating K_(sub) subband signals with output sampling rate f_(sym) and dimension M_(DoF) ×1, where M_(DoF) = M_(sub)M_(feed). In another aspect of the disclosure, the subband generation is performed by passing the ADC output signals provided by each of M_(feed) receiver feeds through a bank of K_(sub) frequency converters and decimators, each with decimated output sampling rate f_(dec) = M_(dec)f_(sym), and then passing each decimated signal through a symbol-rate synchronous channelizer to produce a separate M_(sub) ×1 subband signal with sampling rate f_(sym), on each feed, and then “stacking” those signals to generate a M_(DoF) ×1 vector signal with sampling rate f_(sym), where M_(DoF) = M_(sub)M_(feed). However, the channelizer output data may be generated using many different channelizer structures.

The DSP software can implement algorithms that exploit baseband symbol sequences or navigation data provided by external or “third-party” resources, for each navigation signal of interest to the system, to efficiently obtain RPNT analytics without the need to demodulate the navigation data. Exemplary RPNT analytics can include:

-   estimates of observed received parameters of navigation signals,     including signal geo-observables usable to obtain PNT solutions for     the reception platform, e.g., observed signal frequency-of-arrival     (FOA), time-of-arrival (TOA), and observed local     direction-of-arrival (DOA) or global line-of-bearing (LOB) in     systems employing multi-feed receivers; -   estimates of received signal quality, e.g., received incident power,     and of despread/demodulated navigation sequence quality, e.g.,     despreader output signal-to-interference-and-noise ratio (SINR); and -   measures of accuracy of parameter and signal quality estimates.

RPNT analytic measurement methods include “partially-blind” methods that do not require knowledge of the ranging code for the navigation signals or spatial/polarization array-manifold data (measurement of cross-feed spatial/polarization signatures as a function of DOA) for the receiver, and “copy-aided” methods that may require one or both of the ranging codes or the array manifold to provide more complete RPNT analytics.

In some aspects, referred to herein as “multi-subband” systems, the blind pilot-exploiting adaptive detection and geo-observable estimation DSP algorithm can be implemented over multiple contiguous or non-contiguous channelizer subbands. Within each subband, the method can detect and separate signals without knowledge of the spreading code; in the presence of arbitrary channel multipath, tonal and narrower-band jammers; and in the presence of near-far interference that is much higher than the spreading gain of conventional PNT systems. The method then combines the subbands to provide a detection statistic and geo-observable estimate that exploits the full bandwidth of the system. The method can also provide a natural means for rejecting narrowband jamming, such as by simply “de-emphasizing” subbands containing excessive interference energy. As part of the FOA estimation step, the method may also compute the differential FOA between subbands, thereby allowing direct estimation of the observed transmitter speed. In aspects where multiple navigation signal transmitter ephemeres are available, the method can provide a full PNT solution based on this information alone.

In one aspect, the DSP algorithms can be implemented over multiple portions of the navigation signal band, referred to herein as “navigation signal subbands.” The navigation signal subbands can be obtained using combinations of analog and digital processing. For example, off-center subbands can be obtained by frequency-shifting the SiS with frequency shift f_(R) that differs substantively from the navigation signal transmit frequency, and sampling the frequency-shifted SiS using an ADC with a sampling rate that is much lower than the bandwidth of the navigation signal. Alternately, off-center subbands can be obtained by downconverting, sampling, and frequency-channelizing the navigation signal over its full bandwidth, and then selecting subsets of frequency channels to for each subband. In both cases, the subbands are processed using the adaptive detection, geo-observable estimation, and (for blind methods) symbol stream estimation DSP algorithms described herein. Any combination of methods spanning this range of methods can be similarly employed.

If the number of channelizer DoF’s in each subband is greater than the total number of MOS-DSSS signals (e.g., legitimate navigation signals, spoofers, and substantively time and/or frequency offset multipath images), and tonal interferers in the subband, plus additional processing gain needed to integrate the pilot signal above the non-MOS-DSSS noise-and-interference power level (not limited by the jamming margin of the ranging code), then some aspects of the disclosure can be used to detect and determine geo-observables of all of the MOS-DSSS navigation signals covering that subband without knowledge of the spreading code. Moreover, aspects of the disclosure can be used to perform this detection and geo-observables estimation in the presence of near-far interference that is much higher than the spreading gain of conventional PNT systems, even given the much narrower bandwidth of the single-subband system. In some use scenarios, e.g., reception of navigation signals transmitted from ground or air beacons, where the signal RIP’s are at or near the noise level, little or no additional processing gain is needed to detect the navigation signals, and M_(DoF) can be much less than the spreading gain of the navigation signals. At the same time, the number of baseband symbols needed to obtain a useful detection solution can be greatly reduced by minimizing M_(DoF) . Consequently, this aspect can sharply reduce the complexity of the DSP algorithms, and the susceptibility of the method to “signature blur” induced by motion of the signal transmitters and receivers over the data collection interval.

In some aspects, the navigation signal subbands covered by the channelizer is also dynamically adjusted between collection intervals, e.g., by “stepping” the channelizer through multiple navigation signal subbands, randomly moving the channelizer to different subbands during different collection intervals, or re-tasking the channelizer based on metrics computed during prior collections. This can further reduce susceptibility of the system to narrowband or partial-band jammers attempting to block a portion of the navigation signal band, perhaps a large portion of the band.

In one aspect, dynamic adjustment of the channelizer band is used to identify and remove ambiguities in the FOA estimates provided by the system, based on differential Doppler shift between widely-spaced navigation signal subbands. For example, assuming a 1.2 kilometer/second observed SV speed, the GPS L5 short civil signal can experience a 40 Hz difference in Doppler shift over its ~ 10 MHz bandwidth (80 Hz if the 20 MHz null-to-null bandwidth is assumed), which is well within the resolution capabilities of the method over 1-to-4 second collect intervals, and is sufficient to estimate the observed SV speed, and thereby resolve the 1 kHz FOA ambiguity induced by the 1 ksps baseband symbol rate of this signal. In aspects where multiple navigation signal transmitter ephemeres (time varying position, velocity, acceleration, and transmission time in universal time coordinates) are available, e.g., determined from demodulated navigation data or external resources, methods can provide a full PNT solution based on this information alone.

The DSP detection and geo-observable estimation algorithm can extend without qualitative change to scenarios in which the receiver possesses M_(feed)>1 spatial/polarization diverse receive feeds, except that the channelizer DoF’s are multiplied by the number of receive feeds, such that M_(DoF) = M_(chn) ↺M_(feed), where M_(chn) is the number of channels per feed. However, the resultant system can additionally exploit spatial/polarization diversity of MOS-DSSS signals, including target spoofers that closely emulate the time, frequency, and code diversity of the legitimate navigation signals, and can further excise up to M_(feed)-1 wideband jammers, on the basis of said spatial/polarization diversity, even if those jammers are much stronger than the legitimate navigation signals. As in single-feed aspects, the total DoF’s must additionally be greater than at least the number of MOS-DSSS signals (e.g., legitimate navigation signals and spoofers) and tonal signals impinging on the receiver, over the bandwidth used by the system.

In one aspect of the disclosure, M_(chn) ≈ M_(DoF)/M_(feed), such that M_(DoF) is held constant as the number of feeds is increased. This can further reduce susceptibility of the system to signature blur and provide enhanced rejection of jammers and targeted spoofers attempting to match the time, frequency, and code diversity of legitimate navigation signals, even if those jammers or spoofers are much stronger than the legitimate navigation signals.

In another aspect of the disclosure, a multi-subband channelizer with K_(sub) subbands is implemented on each receiver feed, and M_(chn) and K_(sub) are simultaneously set to M_(chn) ≈ M_(DoF)/M_(feed) and K_(sub) ≈ M_(DoF) ↺M_(feed), such that both M_(DoF) and the total bandwidth used by the system are held constant as the number of feeds is grown. Assuming the constraints on M_(DoF) and M_(feed) given above, the resultant system can achieve the full processing gain of the navigation signal modulation format, and the robustness to near-far interference exhibited by the method in any reception scenario, without unduly affecting complexity or convergence time of the system.

In some aspects of the disclosure, copy-aided post-processing geo-observable estimation methods, along with ranging code information determined from demodulated navigation data or through external resources, are used to resolve the FOA and TOA ambiguities. In systems employing multi-feed receivers, further aspects can use the directional array manifold of the receiver to determine the direction of arrival (DOA) of the received signals; the receiver orientation; and true lines-of-bearing (LOB’s) from the receiver to the signal emitters. In both aspects, the fine geo-observables can be used to geolocate spoofers identified by the system, using known location of the receiver and geo-observables of those spoofers.

In another aspect of the disclosure, the system employs additional network communication methods to allow intercommunication of parameters between multiple receivers. These receivers can employ combinations of single-feed or multi-feed receivers, and can have different ADC rates and channelization methods, as long as each channelizer provides a vector sequence at a 1 ksps output data rate. Among other advantages, this aspect provides additional means for detecting and defeating target spoofers, by using additional network diversity to prevent targeted spoofers from overwhelming any single receiver in the environment.

Any of these aspects can demodulate civil navigation signals in the presence of arbitrary multipath because they do not explicitly exploit the PRN code in the demodulation process. Any of these aspects can provide for determining the yaw-pitch-roll orientation of the receive platform, such as due to relationship between Doppler shift and direction vectors to the GNSS emitters, once the GNSS emitter and receiver positions are determined.

In networked receiver configurations, intercommunication can allow each receiver to continuously draw from and update a “network database” containing estimated TOA’s, FOA’s, and (for multi-feed receivers) DOA’s for navigation signals detected by any receiver in the network, and containing ephemeris information for those receivers. This database can be used to further identify spoofers, including targeted spoofers, based on differential FOA (DFOA) and TOA (DTOA) measured at those receivers, which removes any deliberate timing or frequency shift used by the spoofer to emulate a true navigation signal at a specific victim receiver. If the victim receiver does not employ a multi-feed receiver, and may therefore find itself being successfully spoofed by a targeted attack, this network will provide an unambiguous alert to such an attack, as the spoofed signals will not possess the true GNSS FOA or TOA at the other receivers in the network, and will therefore be removed using the blind adaptive despreading methods employed in aspects of the disclosure.

BRIEF DESCRIPTION OF THE DRAWINGS

Block diagrams and flow diagrams depicted herein can represent computer software instructions or groups of instructions. One or more of the processing blocks or steps may be implemented with modules (e.g., procedures, functions, and so on) that perform the functions described herein. The software codes may be stored in non-transitory computer-readable memory. Alternatively, processing blocks or steps described herein may represent steps performed by hardware, including one or more functionally equivalent circuits, such as a digital signal processor, an application specific integrated circuit, a programmable logic device, a field programmable gate array, a general-purpose processor programmed to perform one or more of the disclosed steps, or other electronic units designed to perform the functions disclosed herein. Disclosed aspects can employ multiple processors communicatively coupled together, and can employ distributed computing, such as Cloud computing, cluster computing, distributed computing, etc.

FIG. 1 depicts a system in which aspects of the disclosure can be implemented.

FIG. 2 illustrates operating parameters and geo-observables that can be employed in aspects disclosed herein.

FIG. 3 is a block diagram that depicts generation of civil navigation signals that aspects of the disclosure can be configured to exploit.

FIG. 4 is a Table of civil GNSS modulation formats that aspects of the disclosure can be configured to exploit.

FIG. 5A, FIG. 5B, FIG. 6A, and FIG. 6B are block diagrams of receiver front-ends that can be configured according to aspects of the disclosure.

FIG. 7 is a block diagram that depicts a single-feed channelizer in accordance with an aspect of the disclosure.

FIG. 8 is a block diagram that depicts a single-feed FFT-based channelizer in accordance with an aspect of the disclosure.

FIG. 9 is a block diagram of a single-feed FFT-based symbol-rate synchronous multi-subband channelizer according to an aspect of the disclosure.

FIG. 10 is a block diagram of a single-feed polyphase filter based symbol-rate synchronous multi-subband channelizer according to an aspect of the disclosure.

FIG. 11 is a block diagram of single-feed symbol-rate synchronous multi-subband channelizer according to another aspect of the disclosure.

FIG. 12A and FIG. 12B depict exemplary channelizer metrics for a rectangularly windowed polyphase filter based symbol-rate synchronous channelizer.

FIG. 13 is a block diagram that depicts a multifeed receiver and channelizer in accordance with an aspect of the disclosure.

FIG. 14 is a block diagram of a multifeed FFT-based symbol-rate synchronous multi-subband channelizer according to an aspect of the disclosure.

FIG. 15 is a block diagram of a multifeed polyphase filter based symbol-rate synchronous multi-subband channelizer according to an aspect of the disclosure.

FIG. 16 and FIG. 17 illustrate despreader parameters that can be employed in aspects disclosed herein.

FIG. 18 is a flow diagram of a multi-subband signal detection, geo-observable/quality estimation, and symbol estimation algorithm that can be employed in aspects disclosed herein.

FIG. 19 is a flow diagram of a multi-subband signal detection and geo-observable/quality estimation, algorithm that can be employed in aspects disclosed herein.

FIG. 20 is a flow diagram of a multi-subband signal detection and geo-observable/quality estimation, algorithm that can be employed in aspects disclosed herein.

FIG. 21 is a flow diagram of a local-maxima search procedure according to an aspect of the disclosure.

FIGS. 22A, 22B, and 22C and FIGS. 23A, 23B, and 23C are plots that depict DFOA search sensitivity.

DETAILED DESCRIPTION

Various aspects of the disclosure are described below. It should be apparent that the teachings herein may be instantiated in a wide variety of forms and that any specific structure, function, or both being disclosed herein are merely representative. Based on the teachings herein one skilled in the art should appreciate that an aspect disclosed herein may be implemented independently of any other aspects and that two or more of these aspects may be combined in various ways. For example, an apparatus may be implemented or a method may be practiced using any number of the aspects set forth herein. In addition, such an apparatus may be implemented or such a method may be practiced using other structure, functionality, or structure and functionality in addition to or other than one or more of the aspects set forth herein.

In the following description, for the purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the invention. It should be understood, however, that the particular aspects shown and described herein are not intended to limit the invention to any particular form, but rather, the invention is to cover all modifications, equivalents, and alternatives falling within the scope of the invention as defined by the claims.

FIG. 1 depicts a model of the presumed normal condition wherein a receiver (1) is receiving GPS signals from an SV (3).

FIG. 2 Illustrates exemplary observable parameters of signals transmitted from a GNSS SV (3) to a GNSS receiver (1) that can be used to geo-locate those signals (“geo-observables”) (6), as a function of parameters local to the GNSS SV (4) and the GNSS receiver (2). These geo-observables (6) include the time-of-arrival (TOA) (after removal of delays induced in transmitter and receiver electronics), received incident power (RIP), line-of-bearing (LOB), direction-of-arrival (DOA) relative to the platform orientation, frequency-of-arrival (FOA) (after removal of known frequency offset between the transmitter and receiver center frequencies), and carrier phase (CCP) (e.g., after calibration and removal of transmitter and receiver phase offsets), which are measurable and differentiable at the receiver (1).

The TOA, RIP, FOA, LOB, DOA, and CCP are referred to as geo-observable, as they are observable parameters of the GNSS received signals that can be used to geo-locate the GNSS receiver, given knowledge of the GNSS transmitter locations over time, which are transmitted within the navigation stream of those signals. They are also referred to here as positioning, navigation, and timing (PNT) analytics, as they provide metadata of the transmitted signals that can be used for purposes beyond geo-location of the receiver, e.g., to assess quality and reliability of the received GNSS signals, in order to aid other PNT systems on board the receiver platform.

FIG. 3 depicts operations used to generate the civil navigation signals that are the focus of the invention. It should be noted that the processing steps shown in FIG. 3 is one conceptual means for generating civil GNSS signals, and ignore other SiS’s that may be transmitted simultaneously with civil GNSS signals, for example, the GPS P, P(Y), M, and L1C SiS’s. In practice, these other SiS’s are largely outside the band occupied by the civil GNSS signals, and moreover are received well below the noise floor. Therefore, their effect on aspects described herein is small. An important exception is military GNSS signals that may be generated by terrestrial pseudolites, as the non-civil components of those signals may be above the noise floor within the civil GNSS passband in some use scenarios.

The processing steps shown in FIG. 3 generate a particular class of spread spectrum signal, referred to herein as a modulation-on-symbol direct-sequence spread spectrum (MOS-DSSS) signal, in which a data symbol sequence d_(T)(n_(sym);ℓ) with symbol rate f_(sym) is spread by a M_(chp) ×1 dimensioned ranging code

c_(T)(𝓁) = [c_(T)(m_(chp); 𝓁)]_(m_(chp) = 0)^(M_(chp) − 1)

that is repeated over every data symbol. The spread signal can be represented as the multiplication of data symbol sequence d_(T)(n_(sym);ℓ) and code vector c_(T)(ℓ) (115), resulting in M_(chp)×1 c_(T)(ℓ)d_(T)(n_(sym);ℓ). The vector signal is then passed through a M_(chp) :1 parallel-to-serial (P/S) converter (116), resulting in a scalar spread signal with rate M_(chp)f_(sym), which can be expressed as c_(T)(n_(chp)modM_(chp);ℓ)d_(T)(n_(chp)/M_(chp);), where (o)modM and denote the modulo-M operation and integer truncation or “floor” operations, respectively. As shall be described below, this property induces inherent massive spectral redundancy to the civil GNSS format that is exploited by disclosed aspects.

The chip-rate spread signal is then passed through digital-to-analog conversion (DAC) (118), and transmission/upconversion stages (120), controlled by an internal clock (clk) (117) and local-oscillator (LO) (119) stages, and transmitted from an antenna (121) to generate a signal-in-space (SiS) given by

$\sqrt{2P_{T}\left( \mathcal{l} \right)}\text{Re}\left\{ {s_{T}\left( {t - \tau_{T}\left( \mathcal{l} \right);\mathcal{l}} \right)e^{j{({2\pi f_{T}t + \varphi_{T}{(\mathcal{l})}})}}} \right\}$

for transmitter ℓ (122), where f_(T) (ℓ) ) is the transmit frequency, typically common to all SV’s (f_(T)(ℓ)≡f_(T)), P_(T)(ℓ), φ_(T)(ℓ) and τ_(T) (ℓ) are the transmit power, phase, and electronics delay for SiS ℓ, respectively, and where s_(T)(t;ℓ) is the complex baseband representation of the transmitted signal.

Based on these operations, the complex-baseband signal can be modeled by

$\begin{matrix} {s_{T}\left( {t;\mathcal{l}} \right) = {\sum\limits_{n_{\text{sym}}}{d_{T}\left( {n_{\text{sym}};\mathcal{l}} \right)h_{T}\left( {t - T_{\text{sym}}n_{\text{sym}};\mathcal{l}} \right)}}\,,} & \text{­­­Eqn (1)} \end{matrix}$

where d_(T) (n_(sym);ℓ) is the encoded symbol stream, T_(sym) = 1/f_(sym), is the symbol period, and h_(T)(t;ℓ) is a spread symbol given by

$\begin{matrix} {h_{T}\left( {t;\mathcal{l}} \right) = {\sum\limits_{m_{\text{chp}} = 0}^{M_{\text{chp}} - 1}{c_{T}\left( {m_{\text{chp}};\mathcal{l}} \right)g_{T}\left( {t - T_{\text{chp}}m_{\text{chp}};\mathcal{l}} \right)}}\,,\quad T_{\text{chp}} = \frac{T_{\text{sym}}}{M_{\text{chp}}},} & \text{­­­Eqn (2)} \end{matrix}$

and where g_(T)(t;ℓ) is the chip symbol modulated by each chip, referred to here as the chip shaping. The spread symbol can also be expressed by it’s continuous Fourier transform H_(T)(f;ℓ) = ∫h_(T)(f;ℓ)e^(-j) ^(2πft)df,

$\begin{matrix} {H_{T}\left( {f;\mathcal{l}} \right) = C_{T}\left( {e^{j2\pi fT_{\text{chp}}};\mathcal{l}} \right)G_{T}\left( {f;\mathcal{l}} \right),} & \text{­­­Eqn (3)} \end{matrix}$

where

$C_{T}\left( {e^{j2\pi f};\mathcal{l}} \right) = {\sum\limits_{m_{\text{chp}}}{c_{T}\left( m_{\text{chp}} \right)e^{- j2\pi fm_{\text{chp}}}}}$

is the discrete Fourier transform of c_(T)(m_(chp);ℓ), and where G_(T)(f;ℓ) is the continuous Fourier transform of g_(T)(t;ℓ). The spread signal can therefore be modeled as a pulse-amplitude modulated (PAM) waveform in which baseband symbol sequence modulates a symbol waveform that is itself direct-sequence spread by the ranging code, leading to the MOS-DSSS nomenclature employed here.

In typical navigation systems, the chip shaping is identical for each transmitted signal (g_(T)(t;ℓ)≡g_(T)(t)); however, in practice, the shaping can also vary between transmitters. For typical navigation systems, the chip shaping is also nominally rectangular, such that G_(T)(f;ℓ)≡T_(chp)Sa(πfT_(chp)), where Sa(x) = sin(x)/x; however, different chip shapings can be employed without changing the basic structure of the MOS-DSSS modulation format.

In typical navigation systems, the spreading code c_(T)(ℓ) employed at transmitter ℓ is unique to that transmitter, and is in fact provisioned by the navigation system. In contrast, the transmit power and center frequency are typically (but not necessarily) identical for each transmitter, while the transmitter carrier phase, delay, and other impairments will vary between the transmitters. However, any of these typical assumptions can be varied in practice, or in different instantiations of a navigation system.

As shown in FIG. 3 , the full symbol sequence d_(T)(n_(sym);ℓ) is generated from an underlying navigation sequence b_(T)(n_(NAV);ℓ)|_(NAV) with navigation rate f_(NAV), which contains information needed to generate a full positioning, navigation, and timing (PNT) solution at navigation receiver. The navigation sequence typically provides timing and ephemeris (time-stamped position, velocity, and acceleration) of the navigation transmitter and larger network, known impairments at the transmitter, and additional information needed to determine quality and timing of the transmitted signal. The navigation sequence also typically provides information about the specific spreading code(s) employed at the transmitter. This navigation sequence is then passed through processing operations to generate an encoded signal sequence b_(T)(n_(sym);ℓ)|_(sym), with symbol rate f_(sym) and for some modulation formats combined with a known pilot sequence p_(T)(n_(sym);ℓ), with symbol rate f_(sym) (112 b) to generate the full symbol sequence signal d_(T)(n_(sym);ℓ). All of these additional processing steps add additional known structure to the transmitted MOS-DSSS signal, which can be exploited in subsequent processing stages. Additionally, the navigation sequence, the information contained within that sequence (e.g., the transmitter ephemeres) can be obtained after transmission from third-party sources, providing additional exploitable structure in applications some aspects of the invention.

FIG. 4 is a tabular description of GNSS networks employing civil GNSS signals that can be modeled as MOS-DSSS signals. As this Figure shows, every GNSS network deployed to date possesses at least one signal type that can be modeled by FIG. 3 . Moreover, all of these signals have a 1 millisecond (ms) ranging code period, and therefore a 1 ksps (kilosample per second) MOS-DSSS baseband symbol rate, albeit with widely varying chip rates and SiS bandwidths, allowing their demodulation using a common receiver structure.

As the last column of FIG. 4 shows, the symbol streams of each of these signals possess properties that can be used to detect and differentiate them from other signals in the environment, including military GNSS signals, jammers, and civil GNSS signals with the same or similar properties. For example, the GPS coarse acquisition (C/A) or “legacy” signal listed in the first row of FIG. 4 is constructed from a 50 bit/second (bps) BPSK (real) navigation sequence b_(T)(n_(NAV);ℓ)|_(NAV), repeated 20 times per MOS-DSSS symbol to generate 1 ksps symbol sequence d_(T)(n_(sym);ℓ) = b_(T)(n_(sym);ℓ)|_(sym) = b_(T)(_(└)n_(sym)/M_(NAV┘);ℓ)|_(NAV) where M_(NAV) = 20. Moreover, b_(T)(n_(NAV):ℓ)|_(NAV) possesses internal fields that can be further exploited to aid the demodulation process, e.g., the TLM Word Preamble transmitted within every 300-bit Navigation subframe.

Similarly, the GPS L5, QZSS, Galileo E5B, E6B, and Beidou 1B transmitters all add a periodic pilot to the quadrature rail of their navigation signals, and the NAVIC RS BOC transmitter adds a periodic pilot to the in-phase rail of its navigation signal. Means for exploiting periodic pilots are a focus of some aspects of this invention.

FIG. 5A and FIG. 5B each depict exemplary direct-conversion front-end reception processing for a single-feed receiver receiving GNSS signals. FIG. 5A shows a direct-conversion receiver front-end that is dedicated to a specific GNSS band, and FIG. 5B shows a direct-conversion receiver front-end that can be flexibly tasked to any GNSS band, consistent with a software defined radio (SDR) architecture. Each can be used as a part of, and in, some of the disclosed aspects.

The direct-conversion receivers shown in FIG. 5A and FIG. 5B each comprise the following:

-   An antenna (128) and low-noise amplifier (LNA) (131), which receives     an SiS containing GNSS emissions within the FoV of the receiver. -   A local oscillator (LO) (133), which generates a complex sinusoid     with frequency f_(R) and phase (φ_(R), comprising an in-phase (I)     real sinusoid, and a quadrature (Q) sinusoid shifted by 90° in phase     from the in-phase sinusoid; -   A complex mixer (134), which multiplies the LNA output signal by the     conjugate of the complex sinusoid, thereby creating a complex signal     with real and imaginary parts, referred to herein as the     in-phase (I) and quadrature (Q) rails of that signal, which has been     downconverted by frequency shift f_(R) -   A dual lowpass filter (LPF) (135), which filters each signal rail to     remove unwanted signal energy from the complex downconverted signal,     yielding analog scalar complex-baseband signal _(xR) (t), centered     at 0 MHz if f_(R) = f_(T) . -   A dual analog-to-digital converter (ADC) (136), which samples each     rail of x_(R)(t) at rate f_(ADC) = M_(ADC)f_(sym), where f_(sym) is     the civil GNSS signal symbol rate shown in FIG. 3 and M_(ADC) is a     positive integer, thereby generating complex dual-ADC output signal     x_(ADC)(n_(ADC)) with time index n_(ADC) (137).

The receiver structure shown in FIG. 5A possesses an additional bandpass filter (BPF) (130), inserted between the antenna and LNA, and dedicated to a specific GNSS band, which removes unwanted adjacent-channel interference (ACI) from the received GNSS SiS, and depicts an LO that is locked to specific, preset value, consistent with a dedicated receiver that is optimized for a specific GNSS band. This structure is more typical of existing GNSS receivers, and provides enhanced signal quality at cost of receiver flexibility.

In contrast, the receiver structure shown in FIG. 5B is consistent with SDR receivers that can be tasked to any GNSS band under user control. This receiver structure can also be used to receive a segment of a GNSS band, e.g., a 1 MHz segment of the GPS L5 band, if f_(R) is significantly offset from the GNSS band center f_(T). The ability to blindly detect, demodulate, and estimate geo-observables of any portion of any GNSS band using the same reception, downconversion, and sampling hardware regardless of the actual bandwidth of the GNSS signals present in that band can be an advantageous feature of at least some of the disclosed aspects.

FIG. 6A and FIG. 6B each depict exemplary front-end reception processing for a single-feed superheterodyne receiver receiving GNSS signals; with FIG. 6A showing a superheterodyne receiver front-end that is dedicated to a specific GNSS band, and FIG. 6B showing a superheterodyne receiver front-end that can be flexibly tasked to multiple GNSS bands. Each can be used as a part of, and in, at least some aspects disclosed herein.

The heterodyne receivers shown in FIG. 6A and FIG. 6B shift the SiS received by an antenna (128) and amplified in an LNA (130) to IF frequency f_(IF), typically using a two-stage mixer 134, 141, yielding analog scalar real-IF signal x_(R)(t) centered on IF frequency f_(IF) if f_(R) = f_(T), with receive phase-shift

φ_(R)=

φ_(R)⁽¹⁾ + φ_(R)⁽²⁾.

The analog signal is then passed to an ADC (143) that samples x_(R)(t) at rate M_(ADC)f_(sym), where f_(sym) is the MOS-DSSS signal symbol rate shown in FIG. 6A and FIG. 6B, and M_(ADC) is a positive integer, yielding real ADC output signal x_(ADC)(n_(ADC)) with time index n_(ADC) (137). The receiver structure shown in FIG. 6A can be optimized for a specific GNSS band, while the receiver structure shown in FIG. 6B can be flexibly tasked to different navigation signal bands, or to segments of such bands, under software control.

In all four cases, the receiver can induce additional delays, frequency offsets, and sampling rate offsets due to internal electronics in the receiver systems, and due to clock rate and timing offset from universal time coordinates (UTC). These are subsumed into the overall channel response, and measured and removed as part of the geo-location solution.

As will be discussed below, the ADC output signal provided by either the direct-conversion or the superheterodyne receivers can be used in disclosed aspects without any change to the design, except control parameters used in channelization operations immediately after that ADC. The disclosed aspects can also be practiced with many other receiver front-ends known to the art. In many cases, the disclosed aspects are blind to substantive differences between those receiver front-ends.

Assuming the reception scenario shown in FIG. 2 , and the direct-conversion receiver structure shown in FIG. 5 , the downconverted and filtered signal shown in x_(R)(t) is given by

$\begin{matrix} {x_{R}(t) = i_{R}(t) + {\sum\limits_{\mathcal{l} = 1}^{L}{s_{R}\left( {t;\mathcal{l}} \right)}}\,,} & \text{­­­Eqn (4)} \end{matrix}$

where i_(R)(t) is the noise and interference received in the receiver passband and s_(R)(t;ℓ) is the complex-baseband representation of the signal received from transmitter ℓ. Over time intervals on the order of 1-10 seconds from reception time t₀, s_(R)(t;ℓ) can be modeled by

$\begin{matrix} \begin{array}{l} {s_{R}\left( {t_{0} + t;\mathcal{l}} \right) \approx g_{TR}\left( \mathcal{l} \right)e^{j2\pi{({\alpha_{TR}{(\mathcal{l})}t + \frac{1}{2}\alpha_{TR}^{(1)}{(\mathcal{l})}t^{2}})}}s_{T}\left( {\left( {1 - \tau_{TR}^{(1)}\left( \mathcal{l} \right)} \right)t -} \right)} \\ {\left( {\tau_{TR}\left( \mathcal{l} \right);\mathcal{l}} \right),} \end{array} & \text{­­­Eqn (5)} \end{matrix}$

at the input to the dual-ADC’s shown in FIG. 5 , where

$\begin{matrix} \begin{array}{l} {g_{TR}\left( \mathcal{l} \right) \triangleq \sqrt{P_{R}\left( {t_{0};\mathcal{l}} \right)}\exp\left( {j\varphi_{TR}\left( \mathcal{l} \right)} \right),\text{end-to-end complex channel}} \\ \text{gain} \end{array} & \text{­­­Eqn (6)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {\varphi_{TR}\left( \mathcal{l} \right) \triangleq \varphi_{T}\left( \mathcal{l} \right) - \varphi_{R} - 2\pi f_{T}\left( \mathcal{l} \right)\tau_{TR}\left( {t_{0};\mathcal{l}} \right),\text{end-to-end channel}} \\ \text{phase} \end{array} & \text{­­­Eqn (7)} \end{matrix}$

$\begin{matrix} {\tau_{TR}\left( \mathcal{l} \right) \triangleq \tau_{T}\left( \mathcal{l} \right) + \tau_{R} + \tau_{TR}^{(0)}\left( {t_{0};\mathcal{l}} \right),\text{end-to-end observed TOA}} & \text{­­­Eqn (8)} \end{matrix}$

$\begin{matrix} {\alpha_{TR}\left( \mathcal{l} \right) \triangleq f_{T}\left( \mathcal{l} \right) + f_{R} - f_{T}\left( \mathcal{l} \right)\tau_{TR}^{(1)}\left( {t_{0};\mathcal{l}} \right),\text{end-to-end observed FOA}} & \text{­­­Eqn (9)} \end{matrix}$

$\begin{matrix} {\tau_{TR}^{(1)}\left( \mathcal{l} \right) \triangleq \tau_{TR}^{(1)}\left( {t_{0};\mathcal{l}} \right),\text{observed differential TOA}\left( \text{DTOA} \right)} & \text{­­­Eqn (10)} \end{matrix}$

$\begin{matrix} {\alpha_{TR}^{(1)}\left( \mathcal{l} \right) \triangleq - f_{T}\left( \mathcal{l} \right)\tau_{TR}^{(2)}\left( {t_{0};\mathcal{l}} \right),\text{observed differential FOA}\left( \text{DFOA} \right).} & \text{­­­Eqn (11)} \end{matrix}$

and where φ_(T)(ℓ), f_(T)(ℓ), and τ_(T)(ℓ) are the carrier phase, frequency, and timing offset induced at transmitter ℓ; φ_(R), f_(R), and τ_(R) are the carrier phase, frequency, and timing offset induced at the receiver; P_(R)(t₀;ℓ) is the received incident power of SiS ℓ at the receiver; and

τ_(TR)^((q))(t₀; 𝓁) = dτ_(TR)(t; 𝓁)/dt^(q)

is the q^(th) derivative of the observed TOA between transmitter ℓ and the receiver at time t. Given observed position, velocity, and acceleration p_(TR)(t;ℓ) = p_(T)(t;ℓ)-p_(R)(t;ℓ), v_(TR)(t;ℓ) = v_(T)(t;ℓ)-v_(R)(t;ℓ), and a_(TR)(t;ℓ) = a_(T)(t;ℓ)-a_(R)(t;ℓ), where

{p_(T)(t; 𝓁), v_(T)(t; 𝓁), a_(T)(t; 𝓁)}_(𝓁 = 1)^(L)

and {p_(R)(t),v_(R)(t),a_(R)(t)} are the transmitter and receiver ephemeres, respectively, then

τ_(TR)^((q))(t; 𝓁)

is given by

$\begin{matrix} {\tau_{TR}^{(0)}\left( {t;\mathcal{l}} \right) = \frac{1}{c}\left\| {\text{p}_{TR}\left( {t;\mathcal{l}} \right)} \right\|_{2},} & \text{­­­Eqn (12)} \end{matrix}$

$\begin{matrix} {\tau_{TR}^{(1)}\left( {t;\mathcal{l}} \right) = \frac{1}{c}\text{u}_{TR}^{T}\left( {t;\mathcal{l}} \right)\text{v}_{TR}\left( {t;\mathcal{l}} \right),} & \text{­­­Eqn (13)} \end{matrix}$

$\begin{matrix} {\tau_{TR}^{(2)}\left( {t;\mathcal{l}} \right) = \frac{1}{c}\left( {\text{u}_{TR}^{T}\left( {t;\mathcal{l}} \right)\text{a}_{TR}\left( {t;\mathcal{l}} \right) + \left\| {\text{P}_{\bot}\left( {\text{p}_{TR}\left( {t;\mathcal{l}} \right)} \right)\text{v}_{TR}\left( {t;\mathcal{l}} \right)} \right\|_{2}} \right),} & \text{­­­Eqn (14)} \end{matrix}$

where u_(TR) (t;ℓ) is the line-of-bearing (LOB) direction vector from the receiver to transmitter ℓ,

$\begin{matrix} {\text{u}_{TR}\left( {t;\mathcal{l}} \right) = \frac{\text{p}_{TR}\left( {t;\mathcal{l}} \right)}{\left\| {\text{p}_{TR}\left( {t;\mathcal{l}} \right)} \right\|_{2}},} & \text{­­­Eqn (15)} \end{matrix}$

and P_(┴) (p_(TR)(t;ℓ)) is the projection matrix orthogonal to p_(TR)(t;ℓ), given by

$\begin{matrix} {\text{P}_{\bot}\left( {\text{p}_{TR}\left( {t;\mathcal{l}} \right)} \right) = \text{I}_{3} - \frac{\text{p}_{TR}\left( {t;\mathcal{l}} \right)\text{p}_{TR}^{T}\left( {t;\mathcal{l}} \right)}{\left\| {\text{p}_{TR}\left( {t;\mathcal{l}} \right)} \right\|_{2}^{2}},} & \text{­­­Eqn (16)} \end{matrix}$

and where

∥^(∘)∥₂

and I₃ denote the L2 Euclidean vector norm and 3 × 3 identity matrix, respectively. As shown in the ‘417 NPA, incorporated herein by reference, the RIP and LOB adheres closely to a zero-order model over time intervals on the order of 1-to-10 seconds, for transmitters in medium-Earth orbits and receivers with dynamics commensurate with high-velocity airborne vehicles.

The local direction-of-arrival (DOA) of the signal received from transmitter ℓ can then be represented by local DOA direction vector u_(R)(t;ℓ) = Ψ_(R)(t)u_(TR)(t;ℓ), which can be converted to azimuth relative to the local receiver heading and elevation relative to the plane of receiver motion using well-known directional transformations. As shown in the ‘417 NPA, incorporated herein by reference, the DOA also adheres closely to a first-order model over time intervals on the order of 1-to-10 seconds, for transmitters in medium-Earth orbits and receivers with dynamics commensurate with high-velocity airborne vehicles.. Moreover, the variation in local DOA is nearly identical for all of the received navigation signals under this scenario, as it primarily due to motion/orientation of the receiver, which equally affects all of the signal LOB’s.

FIG. 7 shows a single-feed receiver structure used in one aspect of the disclosure, describing the symbol-rate synchronous vector channelizer (150) integral to the disclosure. A SiS is received at an antenna (128), and downconverted to complex baseband using a direct-conversion receiver (100) implementing operations such as those shown in FIG. 5A and FIG. 5B, comprising reception, frequency downconversion by receive frequency f_(R) to complex baseband representation, and sampling using a dual ADC (136) with sampling rate M_(ADC)f_(sym), where M_(ADC) is a positive integer. The sampled signal is then passed to a symbol-rate synchronous channelizer (150) comprising the following conceptual operations:

-   A 1:M_(ADC) serial-to-parallel conversion operation that transforms     the rate M_(ADC)f_(sym) ADC output signal to an M_(ADC) ×1 vector     signal x_(sym)(n_(sym)) with rate f_(sym), given by -   $\begin{matrix}     {\text{x}_{\text{sym}}\left( n_{\text{sym}} \right) = \begin{pmatrix}     {x_{R}\left( {T_{\text{sym}}n_{\text{sym}}} \right)} \\      \vdots \\     {x_{R}\left( {T_{\text{sym}}\left( {n_{\text{sym}} + \frac{m_{\text{ADC}} - 1}{M_{\text{ADC}}}} \right)} \right)}     \end{pmatrix}.} & \text{­­­Eqn (17)}     \end{matrix}$ -   An M_(ADC) :K_(chn) channelization operation (155), conceptually     performed by operating on x_(sym)(n_(sym)) using K_(chn) × M_(ADC)     channelization matrix operator T_(chn) (153) to form K_(chn) ×1     channelizer output signal x_(chn)(n_(sym)) = T_(chn)     ∘x_(sym)(n_(sym)) (159), where K_(chn) is the symbol-rate     synchronous output channels used in subsequent DSP operations and     “∘” denotes a general filtering operation. In some aspects of the     invention, T_(chn) (153) is a memoryless matrix and the operation is     a simple matrix multiplication. In other aspect of the invention,     T_(chn) (153) is a more linear-time-invariant (LTI) impulse     response, and the operation is a matrix convolution.

The data rate of x_(chn) (n_(sym)) is nominally equal to the symbol rate of the navigation signal of interest to the system (less offset due to nonideal error in the ADC clock), e.g., 1 ksps for the signals listed in FIG. 4 , and the dimensionality of x_(chn)(n_(sym)) is K_(chn) ×1, where K_(chn) Is the number of frequency channels used by subsequent DSP processing stages. Hence this is referred to here as a symbol-rate synchronous vector channelizer.

A variety of different channelization methods can be implemented using this general processing structure, including time and frequency domain channelization methods, or channelizations employing mixed time-frequency channelizations such as wavelet-based channelizers.

FIG. 8 depicts a specific instantiation of the symbol-rate synchronous vector channelizer (150), using Fast-Fourier Transform (FFT) operations. The channelizer first applies a (nominally zero-padded) discrete Fourier transform (DFT) with M_(ADC) input samples, K_(FFT) ≥ M_(ADC) output bins (152), and a (nominally rectangular) window

{w_(chn)(m_(ADC))}_(m_(ADC) = 0)^(M_(ADC) − 1)

(156) to each S/P output vector, where K_(FFT) can be conveniently chosen to minimize computational complexity and/or cost of the DFT using an FFT algorithm. The channelizer then selects K_(chn) output bins

{k_(bin)(k_(chn))}_(k_(chn) = 0)^(K_(chn) − 1)

for use in subsequent adaptive processing algorithms (154). The channelization operation can be expressed as a multiplication of

$\begin{array}{l} {\text{x}_{\text{sym}}\left( n_{\text{sym}} \right)\text{by}K_{\text{chn}} \times M_{\text{ADC}}\text{matrix T}_{\text{chn}} =} \\ {\left\lbrack {\exp\left( {- j2\pi k_{\text{bin}}\left( k_{\text{chn}} \right)\frac{m_{\text{ADC}}}{K_{\text{FFT}}}} \right)w_{\text{chn}}\left( m_{\text{ADC}} \right)} \right\rbrack} \end{array}$

(153).

In one aspect of the invention, the output channelizer bins are selected contiguously over the passband of the receiver, e.g., using formula k_(bin) (k_(chn)) = (k_(init) +k_(chn))modK_(FFT), where k_(init) is the first FFT bin output from the channelizer (FFT bin associated with channelizer bin 0). As will be discussed herein, an advantage of this approach is minimization of channel signature blur due to movement of the transmitter and/or receiver platforms. In another aspect of the invention, the output channelizer bins are selected sparsely over the passband of the receiver, e.g., using formula k_(bin)(k_(chn)) = (k_(init) +K_(s) _(ep)k_(chn))modK_(FFT), where K_(sep) is an integer separation factor between channelizer bins. An advantage of this approach is ability to avoid known receiver impairments such as LO leakage (e.g., by setting k_(init) to a non-multiple of K_(sep)), or to reduce complexity of the channelization operation, e.g., using sparse FFT methods.

In some aspects of the invention, the output channelizer bins are preselected, e.g., to avoid known receiver impairments such as LO leakage, or to reduce complexity of the channelization operation. In other aspects of the invention, the output channelizer bins are adaptively selected, e.g., to avoid dynamic narrowband co-channel interferers (NBCCI).

In, the ‘417 NPA, which has been incorporated herein by reference, means for processing the vector channelizer output signal, in that case with K_(chn) = M_(DoF)▪ However, exploitation of channelization dimensionalities covering the full active bandwidth of the GPS L1 legacy signal, much less the much wider bandwidth GPS L5 and other civil GNSS signals, would substantively increase both the complexity of the linear-algebraic methods disclosed in the ‘417 NPA, and the reception time needed to implement those methods. For this reason, the ‘417 NPA focused on partial-band and subband methods that only exploited a part of that active bandwidth. Methods for overcoming this limitation without unduly affecting complexity or reception time is a primary focus of this invention.

FIG. 9 depicts an instantiation of a single-feed FFT-based symbol-rate synchronous multi-subband channelizer, intended to overcome these limitations. This instantiation is identical to the FFT-based symbol-rate synchronous channelizer shown in FIG. 8 , except that the K_(chn)-channel selection operation (154) is replaced by a subband selection operation (157) that selects K_(sub) subbands from the K_(FFT) FFT output bins provided by the FFT operation (152). In the aspect of the invention shown in this Figure, the number of bins in each subband is equal to the same number, M_(DoF), referred to herein as the number of degrees of freedom of the subband. This has the benefit of regularizing linear-algebraic DSP operations performed on each subband in subsequent processing steps. However, other aspects of the invention can vary the number of bins, and therefore processing degrees of freedom, e.g., based on numbers of interferers present in each subband or other factors. The FFT output bins selected for the subbands are then organized into K_(sub) M_(DoF) ×1 vector signals given by

$\begin{matrix} \begin{array}{l} {\text{x}_{\text{DoF}}\left( {n_{\text{sym}};k_{\text{sub}}} \right) = \left( \begin{array}{l} {x_{\text{FFT}}\left( {k_{\text{bin}}\left( {0;k_{\text{sub}}} \right),n_{\text{sym}}} \right)} \\  \vdots \\ {x_{\text{FFT}}\left( {k_{\text{bin}}\left( {M_{\text{DoF}} - 1;k_{\text{sub}}} \right),n_{\text{sym}}} \right)} \end{array} \right),k_{\text{sub}} = 0,\ldots,} \\ {K_{\text{sub}} - 1,} \end{array} & \text{­­­Eqn (18)} \end{matrix}$

each with output rate f_(sym), where

{x_(FFT)(k_(bin), n_(sym))}_(k_(bin) = 0)^(K_(FFT) − 1)

are the bins output from the FFT (152).

In one aspect of the invention, the FFT bins selected for each subband are contiguous, e.g., set using formula

k_(bin)(m_(DoF); k_(sub)) = (k_(init)(k_(sub)) + m_(DoF))modK_(FFT).,

where k_(init) (k_(sub)) denotes the initial FFT bin selected for subband k_(sub). As will be discussed herein, an advantage of this approach is minimization of channel signature blur due to movement of the transmitter and/or receiver platforms. In another aspect of the invention, the FFT bins selected for each subband are more sparsely distributed, e.g., using formula

k_(bin)(m_(DoF); k_(sub)) = (k_(init)(k_(sub))) + K_(sep)m_(DoF)modK_(FFT),

where K_(sep) is an integer separation between bins within the subband. An advantage of this approach is ability to avoid known receiver impairments such as LO leakage, or to reduce complexity of the channelization operation, e.g., using sparse FFT methods. The bin spacing K_(sep) can be equal for each subband, or be different between subbands.

In one aspect of the invention, the subband are themselves contiguous, e.g., such that k_(init) (k_(sub)+1) =k_(init) (_(ksub))+M_(DoF))modK_(FFT) for subbands with contiguous bins within each subband, and k_(init) (k_(sub)+1) = (k_(init) (k_(sub))+K_(sep)M_(DoF)) modK_(FFT) for subbands with sparsely separated bins within each subband. In other aspects of the invention, subbands may be widely separated. In yet other aspects of the invention, subbands may be interleaved, e.g., by setting K_(sep)=K_(sub) and k_(init)(k_(sub)) = k_(init)(0)+k_(sub), such that k_(bin)(m_(DoF);k_(sub)) = (k_(init) (0)+(K_(sub)m_(DoF) +k_(sub)))modK_(FFT).

In some aspects of the invention, the output subband bins are preselected, e.g., to avoid known receiver impairments such as LO leakage, or to reduce complexity of the channelization operation. In other aspects of the invention, the output subband bins are adaptively selected, e.g., to avoid dynamic narrowband co-channel interferers (NBCCI).

FIG. 10 is a block diagram of an alternate single-feed polyphase filter based symbol-rate synchronous multi-subband channelizer, according to an aspect of the disclosure. This block diagram also depicts a simple means for separating the channelizer output signal into contiguous subbands. The ADC output data x_(R)(T_(ADC)n_(ADC)) is passed through a 1:M_(ADC) serial-to-parallel converter (151) and a Q_(ADC)M_(ADC)-tap polyphase filter (160) with (typically real, potentially complex) channelizer weights

$\begin{matrix} \begin{array}{l} {\left\{ {w_{\text{chn}}\left( m_{\text{ADC}} \right)} \right\}_{m_{\text{ADC}} = 0}^{Q_{\text{ADC}}M_{\text{ADC}}}(158),\text{given by}} \\ {x_{\text{chn}}\left( {k_{\text{chn}},n_{\text{sym}}} \right)} \\ {= {\sum\limits_{m_{\text{ADC}} = 0}^{Q_{\text{ADC}}M_{\text{ADC}} - 1}{w_{\text{chn}}^{*}\left( m_{\text{ADC}} \right)x_{R}\left( {T_{\text{ADC}}\left( {M_{\text{ADC}}n_{\text{sym}} +} \right)} \right)}}} \\ {\left( \left( m_{\text{ADC}} \right) \right)e^{- j2\pi f_{\text{chn}}{(k_{\text{chn}})}T_{\text{ADC}}m_{\text{ADC}}}} \end{array} & \text{­­­Eqn (19)} \end{matrix}$

such that x_(chn) (n_(sym)) is given by

$\text{x}_{\text{chn}}\left( n_{\text{sym}} \right) = {\sum\limits_{q_{\text{ADC}} = 0}^{Q_{\text{ADC}} - 1}{\text{T}_{\text{chn}}\left( q_{\text{sym}} \right)\text{x}_{\text{sym}}\left( {n_{\text{sym}} + q_{\text{sym}}} \right)}}$

$\begin{matrix} \begin{array}{l} {\text{T}_{\text{chn}}\left( q_{\text{sym}} \right) = \left\lbrack {w_{\text{chn}}^{*}\left( {M_{\text{ADC}}q_{\text{sym}} +} \right)} \right)} \\ {\left( {\left( m_{\text{ADC}} \right)e^{- j2\pi f_{\text{chn}}{(k_{\text{chn}})}T_{\text{ADC}}{({M_{\text{ADC}}q_{\text{sym}} + m_{\text{ADC}}})}}} \right\rbrack.} \end{array} & \text{­­­Eqn (20)} \end{matrix}$

Relating to FIG. 7 , this channelizer employs a transformation x_(chn)(n_(sym)) = T_(chn)◦x_(sym)(n_(sym)), where T_(chn) is a K_(chn)×M_(ADC) polyphase filter with response length Q_(ADC). In the aspect shown in FIG. 10 , the channel frequencies

{f_(chn)(k_(chn))}_(k_(chn) = 0)^(K_(chn) − 1)

are furthermore assumed to satisfy

$\begin{matrix} {f_{\text{chn}}\left( k_{\text{chn}} \right) = \left( {k_{\text{chn}} - \frac{K_{\text{chn}} - 1}{2}} \right)f_{\text{sym}}} & \text{­­­Eqn (21)} \end{matrix}$

such that

$\begin{matrix} \begin{array}{l} {\text{T}_{\text{chn}}\left( q_{\text{sym}} \right) =} \\ {\left( {- 1} \right)^{{({K - 1})}q_{\text{sym}}}\left\lbrack {w_{\text{chn}}^{\ast}\left( {M_{\text{ADC}}q_{\text{sym}} + m_{\text{ADC}}} \right)e^{{- j2\pi k_{\text{chn}}m_{\text{ADC}}}/M_{\text{ADC}}}} \right\rbrack.} \end{array} & \text{­­­Eqn (22)} \end{matrix}$

This channelizer can be implemented using methods well-known to those of ordinary skill in the art.

If K_(chn) = K_(sub)M_(DoF), then x_(chn)(n_(sym)) can be directly transformed to K_(sub) contiguous M_(DoF)×1 vector subband signals

{x_(DoF)(n_(sym); k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1)

with contiguous in-subband channels using a M_(DoF)×K_(sub) reshaping operation (161), such that

x_(DoF)(n_(sym); k_(sub))  =

[x_(chn)(M_(DoF)k_(sub) + m_(DoF), n_(sym))]_(m_(DoF) = 0)^(M_(DoF) − 1).

The subband center frequencies

{f_(sub)(k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1)

are then given by

$\begin{matrix} {f_{\text{sub}}\left( k_{\text{sub}} \right) = \left( {k_{\text{sub}} - \frac{K_{\text{sub}} - 1}{2}} \right)M_{\text{DoF}}f_{\text{sym}},\quad k_{\text{sub}} = 0,\quad,K_{\text{sub}} - 1} & \text{­­­Eqn (23)} \end{matrix}$

while the frequency channels within each subband are locally given by

$\begin{matrix} {f_{\text{chn}}\left( m_{\text{DoF}} \right) = \left( {m_{\text{DoF}} - \frac{M_{\text{DoF}} - 1}{2}} \right)f_{\text{sym}},\quad m_{\text{DoF}} = 0,\quad,M_{\text{DoF}} - 1} & \text{­­­Eqn (24)} \end{matrix}$

such that f_(chn)(M_(DOF)k_(sub) +m_(DoF)) = f_(sub)(k_(sub)) + f_(chn) (m_(DoF)).

FIG. 11 is a block diagram of a single-feed symbol-rate synchronous multi-subband channelizer according to another aspect of the disclosure. In this aspect, the ADC output signal is passed through a set of K_(sub) frequency converters (162 a)-(162 b) and decimators (163 a)-(163 b), each of which frequency downconverts the ADC output signal by frequency f_(sub)(k_(sub)), and decimates that signal to intermediate rate f_(dec) = M_(dec)f_(sym), where M_(dec) is a positive integer. Each decimated signal x_(dec)(n_(dec);k_(sub)) is then passed through a symbol-rate synchronous channelization operation (164 a)-(164 b) comprising a 1:M_(dec) serial-to-parallel converter and a M_(DoF)xM_(dec) channelization operator T_(DoF) (165), to form channelized subband signal x_(DoF)(n_(sym);k_(sub)) = T_(DoF)◦x_(sym)(n_(sym);k_(sub)), where x_(sym) (n_(sym) ;k_(sub)) =

[x_(dec)(M_(dec)n_(sym) + m_(dec); k_(sub))]_(m_(dec) = 0)^(M_(dec) − 1).

The subband frequencies

{f_(sub)(k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1)

can be organized in accordance with Eqn (23); set to a preselected set of arbitrary frequencies; or adaptively determined based on channel dynamics or co-channel interference considerations.

Given the reception scenario described in FIG. 1 and FIG. 2 , and the received signal model given in Eqn (4)-Eqn (14), then the single-sample channelizer output signal shown in FIG. 7 can be modeled by

$\begin{matrix} {\text{x}_{\text{chn}}\left( n_{\text{sym}} \right) = \text{i}_{\text{chn}}\left( n_{\text{sym}} \right) + {\sum\limits_{\mathcal{l} = 1}^{L_{T}}{\text{s}_{\text{chn}}\left( {n_{\text{sym}};\mathcal{l}} \right)}},} & \text{­­­Eqn (25)} \end{matrix}$

where i_(chn)(n_(sym)) = T_(chn)◦i_(sym)(n_(sym)) and s_(chn)(n_(sym);ℓ) = T_(chn)◦s_(sym)(n_(sym);ℓ) are the K_(chn)×1 interference and navigation signal ℓ channelizer output signals, respectively, and where i_(sym()n_(sym)) =

$\begin{array}{l} {\left\lbrack {i_{R}\left( {T_{\text{sym}}\left( {n_{\text{sym}} + \frac{m_{\text{ADC}}}{M_{\text{ADC}}}} \right)} \right)} \right\rbrack\mspace{6mu}\text{and}\mspace{6mu}\text{s}_{\text{sym}}\left( {n_{\text{sym}};\mathcal{l}} \right) =} \\ {\left\lbrack {s_{R}\left( {T_{\text{sym}}\left( {n_{\text{sym}} + \frac{m_{\text{ADC}}}{M_{\text{ADC}}}} \right);\mathcal{l}} \right)} \right\rbrack} \end{array}$

are the M_(ADC)×1 interference and GNSS signal ℓ S/P output vectors, respectively.

The ‘417 NPA, incorporated herein by reference, provides a detailed description of the signals obtaining at the output of a symbol-rate-synchronous channelizer with K_(chn) = M_(DOF). A simpler channel model can be obtained using the polyphase filter aspect shown in FIG. 10 . In this case, under the simplifying assumption

$\begin{matrix} {s_{R}\left( {t;\mathcal{l}} \right) \approx g_{TR}\left( \mathcal{l} \right)e^{j2\pi\alpha_{TR}{(\mathcal{l})}t}s_{T}\left( {t - \tau_{TR}\left( \mathcal{l} \right);\mathcal{l}} \right),} & \text{­­­Eqn (26)} \end{matrix}$

i.e., ignoring affects of platform dynamics other than Doppler shift caused by the transmitter velocity observed at the receiver, then the channelizer output signal s_(chn)(k_(chn),n_(sym);ℓ) is given by

$\begin{matrix} \begin{array}{l} {s_{\text{chn}}\left( {k_{\text{chn}},n_{\text{sym}};\mathcal{l}} \right) =} \\ {e^{j2\pi\alpha_{TR}{(\mathcal{l})}T_{\text{sym}}n_{\text{sym}}}{\sum\limits_{q_{\text{sym}}}{a_{\text{chn}}\left( {k_{\text{chn}},q_{\text{sym}};\mathcal{l}} \right)}}d_{T}\left( {n_{\text{sym}} - q_{\text{sym}};\mathcal{l}} \right),} \end{array} & \text{­­­Eqn (27)} \end{matrix}$

where a_(chn)(k_(chn,)q_(sym);ℓ) is given by

$\begin{matrix} \begin{array}{l} {a_{\text{chn}}\left( {k_{\text{chn}},q_{\text{sym}};\mathcal{l}} \right) = e^{j2\pi{({q_{\text{sym}}T_{\text{sym}} - \tau_{TR}{(\mathcal{l})}})}f_{\text{chn}}{(k_{\text{chn}})}}e^{- j2\pi q_{\text{sym}}\alpha_{TR}{(\mathcal{l})}}} \\ {\times {\int\limits_{- \infty}^{\infty}{W_{\text{chn}}^{\ast}\left( e^{j2\pi fT_{\text{ADC}}} \right)H_{TR}\left( {f + f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)e^{j2\pi{({q_{\text{sym}}T_{\text{sym}} - \tau_{TR}{(\mathcal{l})}})}f}df}},} \end{array} & \text{­­­Eqn (28)} \end{matrix}$

in which w(e^(j2πf)) is the discrete Fourier transform of w_(chn)(m_(ADC)) and H_(TR)(f;ℓ) is the frequency response of the end-to-end symbol shaping, given by

$\begin{matrix} \begin{array}{l} {H_{TR}\left( {f;\mathcal{l}} \right) =} \\ {g_{TR}\left( \mathcal{l} \right)e^{j2\pi\alpha_{TR}{(\mathcal{l})}\tau_{TR}{(\mathcal{l})}}H_{R}(f)H_{T}\left( {f - \alpha_{TR}\left( \mathcal{l} \right);\mathcal{l}} \right)\sqrt{P_{T}\left( \mathcal{l} \right)},} \end{array} & \text{­­­Eqn (29)} \end{matrix}$

where H_(T)(f;ℓ) is the Fourier transform of h_(T)(t;ℓ) given in Eqn (2) and H_(R)(f) is the frequency response of the receiver filtering operations.

If H_(TR)(f;ℓ) is bandlimited below f_(ADC)/2, e.g., using typical pre-ADC antialiasing filters, the bandwidth of W(e^(j) ^(2πf)) is much less than the frequency variation in H_(TR) (f;ℓ), and f_(chn)(k_(chn)) is given by Eqn (21), then Eqn (28) can be approximated by

$\begin{matrix} {a_{\text{chn}}\left( {k_{\text{chn}},q_{\text{sym}};\mathcal{l}} \right) \approx a_{\text{chn}}\left( {k_{\text{chn}};\mathcal{l}} \right)a_{\text{sym}}\left( {q_{\text{sym}};\mathcal{l}} \right),} & \text{­­­Eqn (30)} \end{matrix}$

where a_(chn)(k_(chn);ℓ) and a_(sym)(q_(sym);ℓ) are the frequency-varying and time-varying components of the channel signature,

$\begin{matrix} {a_{\text{chn}}\left( {k_{\text{chn}};\mathcal{l}} \right) = f_{\text{ADC}}H_{TR}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right);\mathcal{l}} \right)e^{- j2\pi f_{\text{chn}}{(k_{\text{chn}})}\tau_{TR}{(\mathcal{l})}},} & \text{­­­Eqn (31)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {a_{\text{sym}}\left( {q_{\text{sym}};\mathcal{l}} \right) =} \\ {\left( {- 1} \right)^{{({K_{\text{chn}} - 1})}q_{\text{sym}}}w_{\text{sym}}^{\ast}\left( {q_{\text{sym}} - \tau_{TR}\left( \mathcal{l} \right)f_{\text{sym}}} \right)e^{- j2\pi q_{\text{sym}}\alpha_{TR}{(\mathcal{l})}},} \end{array} & \text{­­­Eqn (32)} \end{matrix}$

respectively, referred to herein as the channel frequency signature and channel time signature, and where w_(sym)(τ) is the interpolated channelizer window, given by

$\begin{matrix} {w_{\text{sym}}(\tau) = {\sum\limits_{m_{\text{ADC}}}{w_{\text{chn}}\left( m_{\text{ADC}} \right)\text{Sa}\left( {\pi\left( {\tau M_{\text{ADC}} - m_{\text{ADC}}} \right)} \right)}}.} & \text{­­­Eqn (33)} \end{matrix}$

The channelized received signal s_(chn)(k_(chn),n_(sym);ℓ) is then approximated by

$\begin{matrix} {s_{\text{chn}}\left( {k_{\text{chn}},n_{\text{sym}};\mathcal{l}} \right)\quad \approx \quad a_{\text{chn}}\left( {k_{\text{chn}};\mathcal{l}} \right)d_{R}\left( {n_{\text{sym}};\mathcal{l}} \right),} & \text{­­­Eqn (34)} \end{matrix}$

where d_(R)(n_(sym);ℓ) is the signal ℓ symbol sequence observed at the receiver, given by

$\begin{matrix} {d_{R}\left( {n_{\text{sym}};\mathcal{l}} \right)\quad = \quad\left( {a_{\text{sym}}\left( {n_{\text{sym}};\mathcal{l}} \right) \circ d_{T}\left( {n_{\text{sym}};\mathcal{l}} \right)} \right)e^{j2\pi\alpha_{TR}{(\mathcal{l})}T_{\text{sym}}n_{\text{sym}}},} & \text{­­­Eqn (35)} \end{matrix}$

and where “o” here denotes the convolution operation. Further decomposing τ_(TR)(ℓ) and α_(TR)(ℓ) into symbol-normalized TOA and FOA components

$\begin{matrix} \begin{array}{l} {\tau_{TR}\left( \mathcal{l} \right) =} \\ {\left( {{\widetilde{\tau}}_{TR}\left( \mathcal{l} \right) + n_{TR}\left( \mathcal{l} \right)} \right)T_{\text{sym}}\left\{ \begin{array}{ll} {{\widetilde{\tau}}_{TR}\left( \mathcal{l} \right) \in \left\lbrack \begin{array}{ll} 0 & 1 \end{array} \right),} & {\text{­­­Eqn (36)}\mspace{6mu}\text{fine}\mspace{6mu}\text{TOA}} \\ {n_{TR}\left( \mathcal{l} \right) \in {\mathbb{Z}},} & {\text{TOA}\mspace{6mu}\text{symbol}\mspace{6mu}\text{offset},} \end{array} \right)} \end{array} &  \end{matrix}$

$\begin{matrix} \begin{array}{l} {\alpha_{TR}\left( \mathcal{l} \right) =} \\ {\left( {{\widetilde{\alpha}}_{TR}\left( \mathcal{l} \right) + k_{TR}\left( \mathcal{l} \right)} \right)f_{\text{sym}}\left\{ \begin{array}{ll} {{\widetilde{\alpha}}_{TR}\left( \mathcal{l} \right) \in \left\lbrack \begin{array}{ll} {- \frac{1}{2}} & \frac{1}{2} \end{array} \right),} & {\text{­­­Eqn (37)}\mspace{6mu}\text{fine}\mspace{6mu}\text{FOA}} \\ {k_{TR}\left( \mathcal{l} \right) \in {\mathbb{Z}},} & {\text{FOA}\mspace{6mu}\text{Nyquist}\mspace{6mu}\text{zone},} \end{array} \right)} \end{array} &  \end{matrix}$

respectively, then Eqn (34) can be replaced by

$\begin{matrix} {d_{R}\left( {n_{\text{sym}};\mathcal{l}} \right) = \left( {a_{\text{sym}}\left( {n_{\text{sym}};\mathcal{l}} \right) \circ d_{T}\left( {n_{\text{sym}} - n_{TR}\left( \mathcal{l} \right);\mathcal{l}} \right)} \right)e^{j2\pi{\widetilde{\alpha}}_{TR}{(\mathcal{l})}n_{\text{sym}}},} & \text{­­­Eqn (38)} \end{matrix}$

$\begin{matrix} {a_{\text{sym}}\left( {q_{\text{sym}};\mathcal{l}} \right) = \left( {- 1} \right)^{{({K_{\text{chn}} - 1})}q_{\text{sym}}}w_{\text{sym}}^{\ast}\left( {q_{\text{sym}} - {\widetilde{\tau}}_{TR}\left( \mathcal{l} \right)} \right)e^{- j2\pi q_{\text{sym}}\alpha_{TR}{(\mathcal{l})}},} & \text{­­­Eqn (39)} \end{matrix}$

except for phase shift φ(ℓ) = 2π(f_(chn)(0)-ã_(TR)(ℓ))n_(TR)(ℓ), which is subsumed into end-to-end link gain in g_(TR)(ℓ) in Eqn (5). Then d_(R)(n_(sym);ℓ) is completely characterized by the channel time signature defined over fine TOA ranging between 0 and 1; the received TOA measured to integer number of symbols; and the fractional fine FOA ranging between -½ and ½.

For the GNSS signals listed in FIG. 4 , and assuming |a_(TR)(ℓ)| f_(chp), then a_(chn)(k_(chn);ℓ) can be further approximated by

$\begin{matrix} \begin{array}{l} {a_{\text{chn}}\left( {k_{\text{chn}};\mathcal{l}} \right) \approx} \\ {\frac{M_{\text{ADC}}}{\sqrt{M_{\text{chp}}}}\text{Sa}\left( {\pi f_{\text{chn}}\left( k_{\text{chn}} \right)T_{\text{chp}}} \right)H_{R}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)g_{TR}\left( \mathcal{l} \right)\sqrt{P_{T}\left( \mathcal{l} \right)}\varepsilon_{R}\left( {k_{\text{chn}};\mathcal{l}} \right),} \end{array} & \text{­­­Eqn (40)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {\varepsilon_{R}\left( {k_{\text{chn}};\mathcal{l}} \right) =} \\ {\frac{1}{\sqrt{M_{\text{chp}}}}C_{T}\left( {e^{- j2\pi{({f_{\text{chn}}{(k_{\text{chn}})} - \alpha_{TR}{(\mathcal{l})}})}T_{\text{chp}}};\mathcal{l}} \right)e^{- j2\pi{({f_{\text{chn}}{(k_{\text{chn}})} - \alpha_{TR}{(\mathcal{l})};\mathcal{l}})}\tau_{TR}{(\mathcal{l})}},} \end{array} & \text{­­­Eqn (41)} \end{matrix}$

where normalized observed ranging code frequency signature ε_(R)(k_(chn);ℓ) captures the cross-channel frequency variability of each channel frequency signature due to the separate ranging code used at each transmitter, and the differing TOA on each transmission path. For well-designed pseudo-random code sequences and M_(chp) 1, ε_(R)(k_(chn);ℓ) can be modeled as a zero-mean, unit-variance, independent and identically-distributed (i.i.d.) complex-Gaussian random process. The channelized output signal can then be approximated by

$\begin{matrix} {\text{x}_{\text{chn}}\left( n_{\text{sym}} \right) \approx \text{i}_{\text{chn}}\left( n_{\text{sym}} \right) + {\sum\limits_{\mathcal{l} = 1}^{L_{T}}{\text{a}_{\text{chn}}\left( \mathcal{l} \right)d_{R}}}\left( {n_{\text{sym}};\mathcal{l}} \right),} & \text{­­­Eqn (42)} \end{matrix}$

where

a_(chn)(𝓁) = [a_(chn)(k_(chn); 𝓁)]_(k_(chn) = 0)^(K_(chn) − 1)

is the K_(chn)×1 signal ℓ channel frequency signature and a_(chn)(k_(chn);ℓ) is given by Eqn (40)-Eqn (41) and d_(R)(n_(sym);ℓ) is given by Eqn (38).

Assuming the background interference i_(R)(t) given in Eqn (4) is complex-Gaussian and stationary with power spectral density (PSD) S_(iRiR) (f), then the interference power level in each channel is given by

$\begin{matrix} {\left\langle \left| {i_{\text{chn}}\left( {k_{\text{chn}},n_{\text{sym}}} \right)} \right|^{2} \right\rangle = {\int\limits_{{- 1}/2}^{1/2}\left| {W_{\text{chn}}\left( e^{j2\pi f} \right)} \right|^{2}}{\widetilde{S}}_{i_{R}i_{R}}\left( {f + f_{\text{chn}}\left( k_{\text{chn}} \right)T_{\text{ADC}}} \right)df,} & \text{­­­Eqn (43)} \end{matrix}$

$\text{where}\mspace{6mu}{\widetilde{S}}_{i_{R}i_{R}}(f) = f_{\text{ADC}}{\sum\limits_{\mathcal{l}}\left| {H_{R}\left( {\left( {f + \mathcal{l}} \right)f_{\text{ADC}}} \right)} \right|^{2}}S_{i_{R}i_{R}}\left( {\left( {f + \mathcal{l}} \right)f_{\text{ADC}}} \right)$

is the PSD of sampled interference signal i_(R)(T_(ADC)n_(ADC)). Assuming the ADC input signal is bandlimited to less than f_(ADC)/2 ahead of the sampling operation, and that the PSD frequency variation is low over the passband of the polyphase filter frequency response W_(chn)(e^(j2πf)), then Eqn (43) can be approximated by

$\begin{matrix} \begin{array}{l} {\left\langle \left| {i_{\text{chn}}\left( {k_{\text{chn}},n_{\text{sym}}} \right)} \right|^{2} \right\rangle \approx} \\ {f_{\text{ADC}}\left| {H_{R}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)} \right|^{2}S_{i_{R}i_{R}}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right){\int\limits_{{- 1}/2}^{1/2}\left| {W_{\text{chn}}\left( e^{j2\pi f} \right)} \right|^{2}}df} \\ {= f_{\text{ADC}}\left| {H_{R}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)} \right|^{2}S_{i_{R}i_{R}}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)\left\| \text{w}_{\text{chn}} \right\|_{2}^{2}} \\ {= \frac{M_{\text{ADC}}^{2}}{M_{\text{chp}}}f_{\text{chp}}\left| {H_{R}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)} \right|^{2}S_{i_{R}i_{R}}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right),} \end{array} & \text{­­­Eqn (44)} \end{matrix}$

where

$\left\| \text{w}_{\text{chn}} \right\|_{2}^{2} = {\sum\limits_{m_{\text{ADC}}}\left| {w_{\text{chn}}\left( m_{\text{ADC}} \right)} \right|^{2}}$

denotes the squared L2 Euclidean norm of the channelizer weights {w_(chn)(m_(ADC))} (158), and further assuming normalization constraint

∥w_(chn)∥₂² = M_(ADC).

More generally, the interference cross-correlation across symbol lag and frequency channel offset can be approximated by

$\begin{matrix} \begin{array}{l} {R_{i_{\text{chn}}i_{\text{chn}}}\left( {k_{\text{chn}},\mathcal{l}_{\text{chn}};m_{\text{sym}}} \right) \triangleq \left\langle {i_{\text{chn}}\left( {k_{\text{chn}},n_{\text{sym}}} \right)i_{\text{chn}}^{\ast}\left( {\mathcal{l}_{\text{chn}},n_{\text{sym}} - m_{\text{sym}}} \right)} \right\rangle} \\ {\approx \frac{M_{\text{ADC}}^{2}}{M_{\text{chp}}}f_{\text{chp}}S_{i_{R}i_{R}}\left( \frac{f_{\text{chn}}\left( k_{\text{chn}} \right) + f_{\text{chn}}\left( \mathcal{l}_{\text{chn}} \right)}{2} \right)} \\ {\times \left( {- 1} \right)^{{({K_{\text{chn}} - 1})}m_{\text{sym}}}\rho_{\text{chn}}\left( {M_{\text{ADC}}m_{\text{sym}},\frac{k_{\text{chn}} - \mathcal{l}_{\text{chn}}}{M_{\text{ADC}}}} \right),} \end{array} & \text{­­­Eqn (45)} \end{matrix}$

where

$\rho_{\text{chn}}\left( {m,\alpha} \right) = \frac{1}{M_{\text{ADC}}}{\sum\limits_{n}W_{\text{chn}}}(n)W_{\text{chn}}^{*}\left( {n - m} \right)e^{- j2\pi\alpha n}$

is the discrete-time ambiguity function for channelizer window {w_(chn)(m_(ADC))} (158), normalized by M_(ADC) so that |ρ_(chn)(m,α)|≤1.

FIG. 12 depicts the channel time signature (FIG. 12A) and normalized ambiguity function (FIG. 12B) induced by the polyphase-filter based symbol-rate synchronous channelizer (160) shown in FIG. 10 , using single-symbol rectangularly-windowed channelizer weights (158). This Figure also obtains for the FFT-based symbol-rate synchronous channelized shown in FIG. 8 and FIG. 9 . As shown in FIG. 12A, the channel time signature is only substantive over a single symbol for fine TOA offsets between 0 and 1, and is nearly equal to unity within this region. The received symbol sequence d_(R)(n_(sym);ℓ) given in Eqn (38) then simplifies to

$\begin{matrix} {d_{R}\left( {n_{\text{sym}};\mathcal{l}} \right) \approx d_{T}\left( {n_{\text{sym}} - 1 - n_{TR}\left( \mathcal{l} \right);\mathcal{l}} \right)e^{j2\pi{\widetilde{\alpha}}_{TR}{(\mathcal{l})}n_{\text{sym}}},} & \text{­­­Eqn (46)} \end{matrix}$

where the phase term -2πα̃_(TR)(ℓ) is subsumed into g_(TR)(ℓ). Similarly, the normalized ambiguity function shown in FIG. 12B is equal to unity at zero channel and symbol offset. Hence the background interference for this channelizer window can be treated as independent between channels and symbols, yielding K_(chn) × K_(chn) interference autocorrelation matrix (ACM)

$\begin{matrix} \begin{array}{l} {\text{R}_{\text{i}_{\text{chn}}\text{i}_{\text{chn}}} = \frac{M_{\text{ADC}}}{M_{\text{chp}}}f_{\text{chp}}\text{diag}\left\{ {\left| {H_{R}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)} \right|^{2}S_{i_{R}i_{R}}} \right)} \\ {\left( \left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right) \right\}_{k_{\text{chn}} = 0}^{k_{\text{chn}} = - 1}.} \end{array} & \text{­­­Eqn (47)} \end{matrix}$

However, the channel model established above can be extended to any channelizer window.

Multiplying x_(chn)(n_(sym)) by the inverse square-root of R_(ichnichn) yields normalized signal

$\begin{matrix} {{\widetilde{\text{x}}}_{\text{chn}}\left( n_{\text{sym}} \right)\begin{array}{l} {\approx \text{R}_{\text{i}_{\text{chn}}\text{i}_{\text{chn}}}^{- {1/2}}\text{x}_{\text{chn}}\left( n_{\text{sym}} \right),} \\ {\approx \text{ε}_{\text{chn}}\left( n_{\text{sym}} \right) + {\sum\limits_{\mathcal{l} = 1}^{L_{T}}\text{u}_{\text{chn}}}\left( \mathcal{l} \right)d_{R}\left( {n_{\text{sym}};\mathcal{l}} \right),} \end{array}} & \text{­­­Eqn (48)} \end{matrix}$

where d_(R)(n_(sym);ℓ) is given by Eqn (46), ε_(chn)(n_(sym)) is a K_(chn)×1 i.i.d. complex-Gaussian random process with zero mean and ACM I_(Kchn) , and u_(chn)(ℓ) is the K_(chn)×1 signal ℓ normalized channel frequency signature given by

$\begin{matrix} {\text{u}_{\text{chn}}\left( \mathcal{l} \right) = \left\lbrack {\frac{\text{Sa}\left( {\pi f_{\text{chn}}\left( k_{\text{chn}} \right)T_{\text{chp}}} \right)}{\sqrt{f_{\text{chp}}S_{i_{R}i_{R}}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)}}\varepsilon_{R}\left( {k_{\text{chn}};\mathcal{l}} \right)} \right\rbrack g_{TR}\left( \mathcal{l} \right)\sqrt{P_{T}\left( \mathcal{l} \right)}.} & \text{­­­Eqn (49)} \end{matrix}$

In theory, a set of K_(chn)×1 combining weights w̃_(max-SINR)(ℓ) can then be developed that extracts d_(R)(n_(sym);ℓ) from x̃_(chn)(n_(sym)) with maximum-attainable signal-to-interference-and-noise ratio (max-SINR). For the channel model given in Eqn (48)-Eqn (49), and assuming the navigation symbol sequences are independent for each transmitter, w̃_(max-SINR)(ℓ) is given by

$\begin{matrix} {\text{w}_{\text{max-SINR}}\left( \mathcal{l} \right)\begin{array}{l} {\propto \left( {\text{I}_{K_{\text{chn}}} + \text{U}_{\text{chn}}\left( {\sim \mathcal{l}} \right)\text{U}_{\text{chn}}^{H}\left( {\sim \mathcal{l}} \right)} \right)^{- 1}\text{u}_{\text{chn}}\left( \mathcal{l} \right)} \\ \begin{array}{l} {= \left( {\text{I}_{K_{\text{chn}}} - \text{U}_{\text{chn}}\left( {\sim \mathcal{l}} \right)\left( {\text{I}_{L_{T} - 1} + \text{U}_{\text{chn}}^{H}\left( {\sim \mathcal{l}} \right)} \right)^{- 1}\text{U}_{\text{chn}}^{H}} \right)} \\ {\text{u}_{\text{chn}}\left( \mathcal{l} \right),} \end{array} \end{array}} & \text{­­­Eqn (50)} \end{matrix}$

where U_(chn)(~ℓ) = [u_(chn)(ℓ’)]_(ℓ′≠ℓ) is the K_(chn)×(L_(T)-1) matrix of normalized frequency signatures interfering with signal ℓ, and the max-SINR for symbol sequence ℓ is given by

$\begin{matrix} {\gamma_{\text{max-SINR}}\left( \mathcal{l} \right)\begin{array}{l} {\approx \left\| {\text{u}_{\text{chn}}\left( \mathcal{l} \right)} \right\|_{2}^{2}} \\ \begin{array}{l} {- \text{u}_{\text{chn}}^{H}\left( \mathcal{l} \right)\text{U}_{\text{chn}}\left( {\sim \mathcal{l}} \right)\left( {\text{I}_{L_{T} - 1} + \text{U}_{\text{chn}}^{H}\left( {\sim \mathcal{l}} \right)\text{U}_{\text{chn}}\left( {\sim \mathcal{l}} \right)} \right)^{- 1}} \\ {\text{U}_{\text{chn}}^{H}\left( {\sim \mathcal{l}} \right)\text{u}_{\text{chn}}\left( \mathcal{l} \right).} \end{array} \end{array}} & \text{­­­Eqn (51)} \end{matrix}$

Modeling ε_(R)(k_(sub);ℓ) as i.i.d. complex Gaussian with zero mean and unity variance, appropriate for unknown, well-modeled long ranging codes, then Eqn (51) has mean lower bound

$\begin{matrix} {{\overline{\gamma}}_{\text{max-SINR}}\left( \mathcal{l} \right)\begin{array}{l} {= E\left\{ {\gamma_{\text{max-SINR}}\left( \mathcal{l} \right)} \right\}} \\ \begin{array}{l} {\geq \left( {{\overline{\gamma}}_{\text{max-SNR}}\left( \mathcal{l} \right) - \left( {L_{T} - 1} \right)\left( \frac{{\overline{\gamma}}_{\text{max-INR}}\left( \mathcal{l} \right)}{1 + {\overline{\gamma}}_{\text{max-INR}}\left( \mathcal{l} \right)} \right)} \right)} \\ {\left( {\underset{k_{\text{chn}}}{\text{max}}{\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{chn}};\mathcal{l}} \right)} \right),} \end{array} \end{array}} & \text{­­­Eqn (52)} \end{matrix}$

$\begin{matrix} \begin{array}{l} \left. \rightarrow\left\{ \begin{array}{l} {{\overline{\gamma}}_{\text{max-SNR}}\left( \mathcal{l} \right),} \\ {\left( {{\overline{\gamma}}_{\text{max-SNR}}\left( \mathcal{l} \right) - \left( {L_{T} - 1} \right)\max\limits_{k_{\text{chn}}}{\overline{\gamma}}_{\text{max-SNR}}\left( {k_{\text{chn}};\mathcal{l}} \right)} \right),} \end{array} \right) \right. \\ \begin{array}{l} {{\overline{\gamma}}_{\text{max-INR}}\left( \mathcal{l} \right) \ll 1} \\ {{\overline{\gamma}}_{\text{max-INR}}\left( \mathcal{l} \right) \ll 1,} \end{array} \end{array} & \text{­­­Eqn (53)} \end{matrix}$

where γ̅_(Rx-SNR)(k_(chn);ℓ) and γ̅_(max-SNR)(ℓ) are the mean receive SNR (Rx-SNR) and maximum-attainable SNR (max-SNR) of signal ℓ in the absence of interference navigation signals, respectively,

$\begin{matrix} {{\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{chn}};\mathcal{l}} \right)\begin{array}{l} {= E\left\{ \left| {u_{\text{chn}}\left( {k_{\text{chn}};\mathcal{l}} \right)} \right| \right\}} \\ {= \frac{\text{Sa}^{2}\left( {\pi f_{\text{chn}}\left( k_{\text{chn}} \right)T_{\text{chn}}} \right)}{f_{\text{chp}}S_{i_{R}i_{R}}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)}\left| {g_{TR}\left( \mathcal{l} \right)} \right|^{2}P_{T}\left( \mathcal{l} \right),} \end{array}} & \text{­­­Eqn (54)} \end{matrix}$

$\begin{matrix} {{\overline{\gamma}}_{\text{max-SNR}}\left( \mathcal{l} \right)\begin{array}{l} {= E\left\{ \left\| {\text{u}_{\text{chn}}\left( \mathcal{l} \right)} \right\|_{2}^{2} \right\}} \\ {= {\sum\limits_{k_{\text{chn}} = 0}^{K_{\text{chn}} - 1}{{\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{chn}};\mathcal{l}} \right)}},} \end{array}} & \text{­­­Eqn (55)} \end{matrix}$

and where γ̅_(max-INR)(ℓ) is the mean max-attainable SNR of navigation signals interfering with signal ℓ,

$\begin{matrix} {{\overline{\gamma}}_{\text{max-INR}}\left( \mathcal{l} \right) = \frac{1}{L_{T} - 1}{\sum\limits_{\mathcal{l}\prime \neq \mathcal{l}}{{\overline{\gamma}}_{\text{max-SNR}}\left( \mathcal{l}^{\prime} \right).}}} & \text{­­­Eqn (56)} \end{matrix}$

As Eqn (53) shows, at high receive SNR, the max-SINR combiner uses L_(T)-1 combiner degrees of freedom (DoF’s) to excise the navigation signals impinging on the receiver, and employs the remaining K_(chn)-L_(T)+1 combiner DoF’s to suppress the background interference i_(chn)(n_(sym)).

As disclosed in the ‘417 NPA, incorporated herein by reference, linear algebraic methods to determine max-SINR weights can be devised given either knowledge of the content of

{d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)),

e.g., provided by third-party sources; knowledge of a component of

{d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)),

e.g., a known/periodic pilot signal; or knowledge of some other exploitable structure of

{d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)),

as listed in FIG. 4 . These methods can be used to detect the navigation signals in the environment, estimate the max-SINR of the combiner output symbol sequences for those signals, and determine at least the geo-observables

{n_(TR)(𝓁), α̃_(TR)(𝓁)}_(𝓁 = 1)^(L_(T))

for those signals, and if needed estimate information-bearing components of the transmitted symbol streams

{d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)).

Moreover, these techniques can perform these operations without prior knowledge of the signal ranging code or receive frequency signature for the signals, or a search over TOA and FOA is the ranging code is known.

However, as also disclosed in the ‘417 NPA, linear-algebraic methods that can develop combiner weights covering the entire navigation signal passband, e.g., K_(chn) ~ 1,000 for the GPS L1 legacy signal, and K_(chn) ~ 10,000 for the GPS L5 civil signal, require

O(K_(chn)²)

operations/symbol to be implemented. Moreover, they require O(K_(chn)) symbols to converge to a useful solution, e.g., 2-to-4 seconds for the GPS L1 legacy signal, and 20-40 seconds for the GPS L5 civil signal. This convergence time introduces additional complexity in GNSS reception scenarios, due to channel dynamics caused by movement of the GNSS satellite vehicles.

In the ‘417 NPA, partial subband methods that reduce K_(chn) to a manageable value are disclosed to overcome these issues. These methods can provide strong advantage in the presence of strong MOS-DSSS spoofing and jamming signals. However, they sacrifice substantive processing gain to achieve this capability at reasonable complexity and over short reception intervals. This issue can be especially significant for modern wideband navigation signals such as the GPS L5 civil signal. The multi-subband approach provides a path to overcome this limitation.

Returning to FIG. 10 , reshape operation (161) shown in FIG. 10 separates x_(chn)(n_(sym)) into K_(sub) subbands, each comprising MDoF frequency channels such that K_(chn) = K_(sub)M_(DoF), and the signal in each subband can be represented by

$\begin{matrix} {\text{x}_{\text{DoF}}\left( {n_{\text{sym}};k_{\text{sub}}} \right)\begin{array}{l} {= \left( \begin{array}{l} {\text{x}_{\text{chn}}\left( {M_{\text{DoF}}k_{\text{sub}},n_{\text{sym}}} \right)} \\  \vdots \\ {\text{x}_{\text{chn}}\left( {M_{\text{DoF}}k_{\text{sub}} + M_{\text{DoF}} - 1,n_{\text{sym}}} \right)} \end{array} \right)} \\ {\approx \text{i}_{\text{DoF}}\left( {n_{\text{sym}};k_{\text{sub}}} \right) + {\sum\limits_{\mathcal{l} = 1}^{L_{T}}\text{a}_{\text{DoF}}}\left( {k_{\text{sub}};\mathcal{l}} \right)d_{R}\left( {n_{\text{sym}};\mathcal{l}} \right),} \end{array}} & \text{­­­Eqn (57)} \end{matrix}$

where d_(R)(n_(sym;)ℓ), given by Eqn (46) for the rectangularly-windowed channelizer, is shared for all subbands, and a_(DoF)(k_(sub);ℓ) is the M_(DoF)×1signal ℓ frequency signature for subband k_(sub),

$\begin{matrix} {\text{a}_{\text{DoF}}\left( {k_{\text{sub}};\mathcal{l}} \right) = \begin{pmatrix} {a_{\text{chn}}\left( {M_{\text{DoF}}k_{\text{sub}},\mathcal{l}} \right)} \\  \vdots \\ {a_{\text{chn}}\left( {M_{\text{DoF}}k_{\text{sub}} + M_{\text{DoF}} - 1,\mathcal{l}} \right)} \end{pmatrix},} & \text{­­­Eqn (58)} \end{matrix}$

and where i_(DoF)(n_(sym);k_(sub)) is the background noise and interference received in subband k_(sub). If i_(R)(t) is stationary with PSD S_(iR) _(iR) (f) at the LPF output, then the zero-lag ACM of i_(DoF)(n_(sym);k_(sub)) is further approximated by

$\begin{matrix} \begin{array}{l} {\text{R}_{\text{i}_{\text{DoF}}\text{i}_{\text{DoF}}}\left( k_{\text{sub}} \right) \approx \frac{M_{\text{ADC}}^{2}}{M_{\text{chp}}}f_{\text{chp}}\left| {H_{R}\left( {f_{\text{sub}}\left( k_{\text{sub}} \right)} \right)} \right|^{2}} \\ {S_{i_{R}i_{R}}\left( {f_{\text{sub}}\left( k_{\text{sub}} \right)} \right)\text{P}_{\text{sub}},} \end{array} & \text{­­­Eqn (59)} \end{matrix}$

where

$f_{\text{sub}}\left( k_{\text{sub}} \right) = \left( {k_{\text{sub}} - \frac{K_{\text{sub}} - 1}{2}} \right)M_{\text{DoF}}f_{\text{sym}}$

is the center frequency of subband k_(sub′) and where

P_(sub)

$= \left\lbrack {\frac{1}{M_{\text{ADC}}}{\sum\limits_{m}\left| {w_{\text{chn}}(m)} \right|^{2}}e^{- j2\pi{({k - \mathcal{l}})}{m/M_{\text{ADC}}}}} \right\rbrack_{k,\mathcal{l} = 0}^{M_{\text{DoF}} - 1}$

is the normalized subband interference ACM, which is equal to I_(MDoF) for the single-symbol rectangular channelizer window.

If M_(DoF) ≥ L_(T), then navigation symbol sequences ℓ can be received in each subband using an appropriately designed set of max-SINR linear in-subband linear combining weights w_(max-SINR)(k_(sub);ℓ), such that combiner output symbol sequence

d̂_(R)(n_(sym); k_(sub); 𝓁) = w_(max-SINR)^(H)(k_(sub); 𝓁)x_(sub)(n_(sym); k_(sub))

excises the interfering navigation signals {d_(R)(n_(sym);ℓ’)}_(ℓ′≠ℓ) also received in the subband using L_(T)-1 degrees of freedom, and suppresses the background interference i_(sub)(n_(sym);k_(chn)) in the subband using the remaining M_(DoF)-L_(T)+1 degrees of freedom, and where (↺)^(H) denotes the Hermitian transpose operation. Assuming a nearly-flat signature and interference response over each subband, and rectangular channelizer windows, then the SINR achievable using this linear combiner satisfies mean lower bound

$\begin{matrix} \begin{array}{l} {{\overline{\gamma}}_{\text{max-SINR}}\left( {k_{\text{sub}};\mathcal{l}} \right) \geq \left( {M_{\text{DoF}} - \left( {L_{T} - 1} \right)\left( \frac{{\overline{\gamma}}_{\text{max-INR}}\left( {k_{\text{sub}};\mathcal{l}} \right)}{1 + {\overline{\gamma}}_{\text{max-INR}}\left( {k_{\text{sub}};\mathcal{l}} \right)} \right)} \right)} \\ {{\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{sub}};\mathcal{l}} \right),} \end{array} & \text{­­­Eqn (60)} \end{matrix}$

$\begin{matrix} \left. \rightarrow\left\{ {\begin{array}{l} {M_{\text{DoF}}{\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{sub}};\mathcal{l}} \right),} \\ {\left( {\text{M}_{\text{DoF}} - L_{T} + 1} \right){\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{sub}};\mathcal{l}} \right),} \end{array}\begin{array}{l} {{\overline{\gamma}}_{\text{max-INR}}\left( {k_{\text{sub}};\mathcal{l}} \right) \ll 1} \\ {{\overline{\gamma}}_{\text{max-INR}}\left( {k_{\text{sub}};\mathcal{l}} \right) \gg 1} \end{array}} \right) \right. & \text{­­­Eqn (61)} \end{matrix}$

in each subband, where γ̅_(Rx-SNR)(k_(sub);ℓ) and γ̅_(max-INR)(k_(sub);ℓ) are the mean receive signal ℓ SNR and the max-attainable navigation signal interference received along with that signal in subband k_(sub′) respectively

$\begin{matrix} {{\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{sub}};\mathcal{l}} \right) \approx \frac{\text{Sa}^{2}\left( {\pi f_{\text{sub}}\left( k_{\text{sub}} \right)T_{\text{chp}}} \right)}{f_{\text{chp}}S_{i_{R}i_{R}}\left( {f_{\text{sub}}\left( k_{\text{sub}} \right)} \right)}\left| {g_{TR}\left( \mathcal{l} \right)} \right|^{2}P_{T}\left( \mathcal{l} \right),} & \text{­­­Eqn (62)} \end{matrix}$

$\begin{matrix} {{\overline{\gamma}}_{\text{max-INR}}\left( {k_{\text{sub}};\mathcal{l}} \right) \approx \frac{M_{\text{DoF}}}{L_{T} - 1}{\sum\limits_{\mathcal{l}^{\prime} \neq \mathcal{l}}{\overline{\gamma}}_{\text{Rx-SNR}}}\left( {k_{\text{sub}};\mathcal{l}^{\prime}} \right).} & \text{­­­Eqn (63)} \end{matrix}$

As Eqn (61) shows, the max-SINR weights can excise up to L_(T)-1 strong MOS-DSSS signals impinging on the receiver in each subband, and can provide an additional factor of M_(DoF)-L_(T)+1 SNR gain to suppress the remaining background noise and interference in each subband. In low-SNR environments where all of the navigation signals are received well below the noise floor (γ̅_(Rx-SNR)(k_(sub);ℓ)«⅟M_(DoF)), the max-SINR weights can provide the full factor-of-M_(DoF) SNR gain available to the receiver.

In aspects of the invention, the symbol sequence estimates from each subband are then further combined using cross-band linear combining weights

{g_(sub)(k_(sub); 𝓁)}_(K_(sub) = 0)^(K_(sub) − 1),

yielding multi-subband symbol sequence estimate

$\begin{matrix} {{\hat{d}}_{R}\left( {n_{\text{sym}};\mathcal{l}} \right)\begin{array}{l} {\approx {\sum\limits_{k_{\text{sub}} = 0}^{K_{\text{sub}} - 1}g_{\text{sub}}^{\ast}}\left( {k_{\text{sub}};\mathcal{l}} \right){\hat{d}}_{R}\left( {n_{\text{sym}};k_{\text{sub}};\mathcal{l}} \right)} \\ {\approx \text{g}_{\text{sub}}^{H}\left( \mathcal{l} \right){\hat{\text{d}}}_{\text{sub}}\left( {n_{\text{sym}};\mathcal{l}} \right).} \end{array}} & \text{­­­Eqn (64)} \end{matrix}$

Assuming each subband has substantively excised the MOS-DSSS signals impinging on the receiver, and has provided a unit-power output signal, then d _(sub)(n_(sym);ℓ) can be approximated by

${\hat{\text{d}}}_{R}\left( {n_{\text{sym}};\mathcal{l}} \right) \approx \Sigma_{\text{max-SINR}}^{1/2}\left( \mathcal{l} \right)\varepsilon\left( {n_{\text{sym}};\mathcal{l}} \right) + \text{a}{}_{\text{max-SINR}}\left( \mathcal{l} \right){\hat{d}}_{R}\left( {n_{\text{sym}};\mathcal{l}} \right)$

$\Sigma_{\text{max-SINR}}\left( \mathcal{l} \right) = \text{diag}\left\{ \frac{\gamma_{\text{max-SINR}}\left( {k_{\text{sub}};\mathcal{l}} \right)}{\left( {1 + \gamma_{\text{max-SINR}}\left( {k_{\text{sub}};\mathcal{l}} \right)} \right)^{2}} \right\}$

$\text{a}_{\text{max-SINR}}\left( \mathcal{l} \right) = \left\lbrack \frac{\gamma_{\text{max-SINR}}\left( {k_{\text{sub}};\mathcal{l}} \right)}{1 + \gamma_{\text{max-SINR}}\left( {k_{\text{sub}};\mathcal{l}} \right)} \right\rbrack.$

The SINR of this signal is maximized by setting g_(sub)(ℓ) = [1+γ_(Max-SINR)(k_(sub);ℓ)], yielding output SINR

$\begin{matrix} {\gamma_{\text{max-SINR}}\left( {n_{\text{sym}};\mathcal{l}} \right) = {\sum\limits_{k_{\text{sub}} = 0}^{K_{\text{sub}} - 1}{\gamma_{\text{max-SINR}}\left( {k_{\text{sub}};\mathcal{l}} \right),}}} & \text{­­­Eqn (65)} \end{matrix}$

which has mean lower bound

$\begin{matrix} \begin{array}{l} {{\overline{\gamma}}_{\text{max-SINR}}\left( \mathcal{l} \right) \geq {\sum\limits_{k_{\text{sub}} = 0}^{K_{\text{sub}} - 1}\left( {M_{\text{DoF}} - \left( {L_{T} -} \right)} \right)}} \\ {\left( {(1)\left( \frac{{\overline{\gamma}}_{\text{max-INR}}\left( {k_{\text{sub}};\mathcal{l}} \right)}{1 + {\overline{\gamma}}_{\text{max-INR}}\left( {k_{\text{sub}};\mathcal{l}} \right)} \right)} \right){\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{sub}};\mathcal{l}} \right),} \end{array} & \text{­­­Eqn (66)} \end{matrix}$

$\begin{matrix} \left. \rightarrow\left\{ \begin{array}{ll} {M_{\text{DoF}}{\sum\limits_{k_{\text{sub}}}{{\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{sub}};\mathcal{l}} \right)}},} & {{\overline{\gamma}}_{\text{max-INR}}\left( {k_{\text{sub}};\mathcal{l}} \right) \ll 1} \\ {\left( {M_{\text{DoF}} - L_{T} + 1} \right){\sum\limits_{k_{\text{sub}}}{{\overline{\gamma}}_{\text{Rx-SNR}}\left( {k_{\text{sub}};\mathcal{l}} \right)}},} & {{\overline{\gamma}}_{\text{max-INR}}\left( {k_{\text{sub}};\mathcal{l}} \right) \gg 1.} \end{array} \right) \right. & \text{­­­Eqn (67)} \end{matrix}$

Thus, the multi-subband solution can recover all or most of the processing gain of the full-channelizer max-SINR combiner, if M_(DoF) is greater than the expected number of MOS-DSSS signals impinging on the receiver.

In-subband linear combining weights approaching the max-SINR solution can be computed using linear-algebraic methods disclosed in the ‘417 NPA, given either knowledge of the content of

{d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)),

e.g., provided by third-party sources; knowledge of a component of

{d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)),

e.g., a known/periodic pilot signal; or other exploitable structure of

{d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)),

as listed in FIG. 4 . These methods also can be used to estimate the max-SINR obtaining in each subband, thereby providing cross-subband weights for the second combining stage shown in Eqn (64). The multi-subband solution can further be used to estimate the full max-SINR of the received symbol sequences, detect the navigation signals in the environment, determine at least the geo-observables

{n_(TR)(𝓁), α̃_(TR)(𝓁)}_(𝓁 = 1)^(L_(T))

for those signals, and if needed estimate information-bearing components of the transmitted symbol streams

{d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)).

Moreover, these techniques can perform these operations without prior knowledge of the signal ranging code or receive frequency signature for the signals, or a search over TOA and FOA is the ranging code is known.

In contrast to the full-band solution, however, the complexity of the multi-subband method scales quadatrically with MDoF and linearly with K_(sub), such that the full set of K_(chn) channels provided by the symbol-rate synchronous channelizer can be processing at a complexity of

O(K_(sub)M_(DoF)²) = O(K_(chn)M_(DoF))

operations/symbol – a factor of K_(sub) saving in operations relative to the full-band system. Moreover, the number of symbols needed to provide a stable solution using the methods disclosed in the ‘417 NPA are O(M_(DoF)), also a factor of K_(sub) reduction in convergence time relative to the full-band methods. In the invention, the value of M_(DoF) becomes a design parameter dictated by the number of MOS-DSSS signals expected to be encountered in the receive environment, and the desired complexity and convergence time of the signal reception algorithms, independent from the bandwidth of the receiver.

As an specific example, setting K_(chn) =960 for the GPS L1 legacy signal and K_(chn) =9,600 for the GPS L5 civil signal, i.e., and setting M_(DoF) =60 in both cases to allow reception of all of the GPS signals in the field of view of the receiver as well as 30-45 beacons or spoofers, then the multi-subband processor has a complexity and convergence time reduction of 16 for the GPS L1 legacy signal, and 160 for the GPS L5 civil signal, allowing convergence in as little as 120 ms for both methods and complexity of O(57.6) Mops/s for the GPS L1 legacy signal and O(576) Mops/s for the GPS L5 civil signal – well within the capabilities of modern DSP equipment, and low-cost DSP for the GPS L1 legacy signal. Moreover, the multi-subband approach is highly amenable to parallel processing methods, facilitating it’s implementation in FPGA’s and general-purpose GPU (GPGPU) architectures.

Taking channel time-variation given in Eqn (5)-Eqn (14) due to moving transmitters and/or receivers more precisely into account, then over reception intervals on the order of 0-to-2 seconds, Eqn (42) generalizes to

$\begin{matrix} {\text{x}_{\text{chn}}\left( n_{\text{sym}} \right) \approx \text{i}_{\text{chn}}\left( n_{\text{sym}} \right) + {\sum\limits_{\mathcal{l} = 1}^{L_{\text{T}}}{\text{a}_{\text{chn}}\left( {n_{\text{sym}};\mathcal{l}} \right)d_{R}\left( {n_{\text{sym}};\mathcal{l}} \right)}}\,,} & \text{­­­Eqn (68)} \end{matrix}$

where time-varying frequency-signature a_(chn)(n_(sym);ℓ) is approximated by

$\begin{matrix} {\text{a}_{\text{chn}}\left( {n_{\text{sym}};\mathcal{l}} \right) \approx \text{a}_{\text{chn}}\left( \mathcal{l} \right) \odot \left\lbrack e^{- j2\pi{({f_{\text{chn}}{(k_{\text{chn}})} - \alpha_{TR}{(\mathcal{l})}})}T_{\text{sym}}{({{\widetilde{\tau}}_{TR}{(\mathcal{l})} + \tau_{TR}^{(1)}{(\mathcal{l})}n_{\text{sym}}})}} \right\rbrack,} & \text{­­­Eqn (69)} \end{matrix}$

and the observed symbol sequence is given by

$\begin{matrix} \begin{array}{l} {d_{R}\left( {n_{\text{sym}};\mathcal{l}} \right) \approx \left( {a_{\text{sym}}\left( {n_{\text{sym}};\mathcal{l}} \right) \otimes d_{T}\left( {n_{\text{sym}} -} \right)} \right)} \\ {\left( \left( {n_{TR}\left( \mathcal{l} \right);\mathcal{l}} \right) \right)e^{j2\pi{({{\widetilde{\alpha}}_{TR}{(\mathcal{l})}n_{\text{sym}} + \frac{1}{2}{\widetilde{\alpha}}_{TR}^{(1)}{(\mathcal{l})}n_{\text{sym}}^{2}})}}} \\ \left. \rightarrow d_{T}\left( {n_{\text{sym}} - n_{TR}\left( \mathcal{l} \right) - 1;\mathcal{l}} \right)e^{j2\pi{({{\widetilde{\alpha}}_{TR}{(\mathcal{l})}n_{\text{sym}} + \frac{1}{2}{\widetilde{\alpha}}_{TR}^{(1)}{(\mathcal{l})}n_{\text{sym}}^{2}})}}, \right. \end{array} & \text{­­­Eqn (70)} \end{matrix}$

and where “ ” denotes the element-wise multiplication (Hadamard product) operation,

τ_(TR)⁽¹⁾(𝓁)

is the observed differential TOA (DTOA) given in Eqn (10), and

α̃_(TR)⁽¹⁾(𝓁) = α_(TR)⁽¹⁾(𝓁)T_(sym)²

is the symbol-normalized observed differential FOA (DFOA). Hence the DTOA principally affects the channel frequency signature, by inducing a linear frequency shift across the symbol-rate synchronous frequency channels, while the DFOA principally affects the channel time signature, by inducing a quadratic frequency shift over the time symbols. However, it should be noted that the FOA is linearly related to the DTOA, as shown in Eqn (9), hence it affects both signature components.

For typical GNSS reception scenarios,

τ_(TR)⁽¹⁾(𝓁)

is typically less than 4 µs/s in magnitude, resulting in a link FOA of ±6.3 kHz at the L1 GPS center frequency, and ±4.7 kHz at the GPS L5 center frequency. From Eqn (69), This DTOA value can also cause a differential frequency shift of as much as 8 Hz over the ±1.023 MHz L1 GPS legacy signal band, or 80 Hz over the ±10.23 MHz L5 civil GPS band; hence it must be taken into account in any geo-observable estimation algorithm implemented over > 10 ms reception intervals.

In this regard, the multi-subband approach shown in FIG. 10 provides an additional benefit, by separating the frequency channels into narrow subbands during the reshape operation (161). In this case, the subband signal can be modeled by

$\begin{matrix} {\text{x}_{\text{DoF}}\left( n_{\text{sym}} \right) \approx \text{i}_{\text{DoF}}\left( n_{\text{sym}} \right) + {\sum\limits_{\mathcal{l} = 1}^{L_{\text{T}}}{\text{a}_{\text{DoF}}\left( {n_{\text{sym}};\mathcal{l}} \right)d_{R}\left( {n_{\text{sym}};k_{\text{sub}};\mathcal{l}} \right)}}\,,} & \text{­­­Eqn (71)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {d_{R}\left( {n_{\text{sym}};k_{\text{sub}};\mathcal{l}} \right) = d_{T}\left( {n_{\text{sym}} - n_{TR}\left( \mathcal{l} \right) -} \right)} \\ {\left( {1;\mathcal{l}} \right)e^{- j2\pi{({f_{T} + f_{\text{sub}}{(k_{\text{sub}})}})}T_{\text{sym}}{({\tau_{TR}^{(1)}{(\mathcal{l})}n_{\text{sym}} + \frac{1}{2}{\widetilde{\tau}}_{TR}^{(2)}{(\mathcal{l})}n_{\text{sym}}^{2}})}},} \end{array} & \text{­­­Eqn (72)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {\text{a}_{\text{DoF}}\left( {n_{\text{sym}};\mathcal{l}} \right) \approx \text{a}_{\text{DoF}}\left( \mathcal{l} \right) \odot} \\ {\left\lbrack e^{- j2\pi{({m_{\text{DoF}} - {(\frac{M_{\text{DoF}} - 1}{2})} - \alpha_{TR}{(\mathcal{l})}T_{\text{sym}}})}{({{\widetilde{\tau}}_{TR}{(\mathcal{l})} + \tau_{TR}^{(1)}{(\mathcal{l})}n_{\text{sym}}})}} \right\rbrack.} \end{array} & \text{­­­Eqn (73)} \end{matrix}$

In aspects of the invention described herein, the substantive cross-subband frequency shift modeled in Eqn (72) is estimated as part of the geo-observable estimation procedure. The observed time-variability in the subband frequency signature a_(DoF)(n_(sym);ℓ) modeled in Eqn (73) induces a much lower in-subband frequency shift, due to the much lower bandwidth of the subband. For the example cited above, where M_(DoF) = 60, a DTOA of 4 µs/s induces a frequency shift of 0.24 Hz across the subband, or 0.12 cycles over 500 ms. As disclosed in the ‘417 NPA, this frequency shift creates low-level dispersive components that can be excised as part of the adaptation process for sufficiently large values of M_(DoF); the value considered here is likely to be more than sufficient for typical GNSS reception scenarios, even in the presence of strong ground beacons, and in the presence of strong spoofers attempting to emulate GNSS transmitter dynamics.

Similarly, for typical GNSS reception scenarios,

α_(TR)⁽¹⁾(𝓁)

is typically less than 3 Hz/s in magnitude. As a consequence, its affect is minimal for reception intervals of 500 ms or less, and manageable for reception intervals on the order of 2 seconds or less, using aspects of the invention disclosed herein.

Accounting for worst-case channel dispersion, if M_(DoF) ≥ 2L_(T) the entire network symbol stream can be extracted from the channel, i.e., the interfering MOS-DSSS symbol sequences can be excised from each GNSS symbol, using purely linear combining operations, and even if the navigation signals are received with high-SNR inducing significant cross-interference. This can include methods well known to the signal processing community, e.g., blind adaptive baseband extraction methods described in the prior art, and linear minimum-mean-square-error (LMMSE) methods described in the prior art. For signals transmitted from MEO GNSS SV’s, and in the absence of pseudolites or spoofers operating in the same band as those SV’s, no more than 12 SV’s are likely to be within the field of view of a GNSS receiver at any one time (L_(T) ≤ 12). Over short observation intervals, maximum substantive rank change induced by the MOS-DSSS signals is therefore likely to be 24. This number can grow to 48 in the presence of 12 spoofers “assigned” to each legitimate navigation signal. In both cases, this is much less than the number of channels available for any GNSS signals listed in FIG. 4 .

This channel response is also exactly analogous to channel responses induced in massive MIMO networks currently under investigation for next-generation (5G) cellular communication systems. However, it achieves this response using only a single-feed receiver front-end, thereby bypassing the most challenging aspect of massive MIMO transceiver technology. And it provides an output signal with an effective data rate of 1 ksps for all of the signals listed in FIG. 4 , as opposed to 10-20 Msps for massive MIMO networks under consideration for 5G cellular communication applications. This receiver front-end also provides considerable flexibility in the number of degrees of freedom M_(DoF) used at the receiver, e.g., by modifying the ADC sampling rate, or the number of FFT bins used in the channelizer, or any combination thereof (albeit, at loss in processing gain if the channelizer does not fully cover the GNSS signal passband).

Lastly, the digital signal processing (DSP) operations needed to exploit this channel response are expected to be very similar to operations needed for 5G data reception, albeit at a 3-4 order-of-magnitude lower switching rate. Given the massive investment expected in 5G communications over the next decade, and the ongoing exponential improvements in cost and performance of DSP processing and memory, e.g., Moore’s and Kryder’s Laws, the ability to fully exploit this channel response will become increasingly easier over time.

The previous description is provided to enable any person skilled in the art to practice the various aspects described herein. Various modifications to these aspects will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other aspects. Thus, the claims are not intended to be limited to the aspects shown herein, but is to be accorded the full scope consistent with the language claims, wherein reference to an element in the singular is not intended to mean “one and only one” unless specifically so stated, but rather “one or more.” The word “exemplary” is used herein to mean “serving as an example, instance, or illustration.” Any aspect described herein as “exemplary” is not necessarily to be construed as preferred or advantageous over other aspects. Unless specifically stated otherwise, the term “some” refers to one or more. The phrase “A or B” may correspond to A only, B only, or A and B. Combinations such as “at least one of A, B, or C,” “one or more of A, B, or C,” “at least one of A, B, and C,” “one or more of A, B, and C,” and “A, B, C, or any combination thereof” include any combination of A, B, and/or C, and may include multiples of A, multiples of B, or multiples of C. Specifically, combinations such as “at least one of A, B, or C,” “one or more of A, B, or C,” “at least one of A, B, and C,” “one or more of A, B, and C,” and “A, B, C, or any combination thereof” may be A only, B only, C only, A and B, A and C, B and C, or A and B and C, where any such combinations may contain one or more member or members of A, B, or C. The words “module,” “block,” “element,” “device,” and the like may not be a substitute for the word “means.” As such, no claim element is to be construed as a means plus function unless the element is expressly recited using the phrase “means for.”

EXTENSION TO MULTIFEED RECEPTION

FIG. 13 depicts aspects of the invention, in which the MOS-DSSS navigation signals are received over M_(feed) feeds using a set of spatial and/or polarization diverse antennas (90 a-90 c). In some aspects, that antennas are additionally coupled from a larger set of antennas to M_(feed) feeds using a beamforming network (BFN) (94), e.g., a Butler mode former, to provide spatial beamsteering operations, exclude known interference sources, or extract stronger and/or spatially whitened modes from the antenna array. The M_(feed) feeds are then coherently converted to complex-baseband using a direct-conversion receiver bank (200), to form

M_(feed) × 1signal vector x_(R)(t) = [x_(R)(t; m_(feed))]_(m_(feed) = 1)^(M_(feed)).

The feed vector x_(R)(t) is then sampled by a bank of dual-ADC’s (243), at sampling rate M_(ADC)f_(sym) determined by a common clock (117) where f_(sym) is the baseband symbol rate of the MOS-DSSS navigation signals transmitted to the receiver, and M_(ADC) is a positive integer.

Assuming coherent downconversion of the received signals, the complex M_(feed) ×1 signal generated output from the receiver bank (200) is modeled as

$\text{x}_{R}(t) = \text{i}_{R}(t) + {\sum\limits_{\mathcal{l} = 1}^{L}{\text{s}_{R}\left( {t;\mathcal{l}} \right)}},$

where i_(R)(t) comprises the M_(feed) ×1 vector of background noise and co-channel interference (CCI) present in the receiver passband, and S_(R)(t;ℓ) is M_(feed) ×1 vector of MOS-DSSS navigation signals received from transmitted ℓ. In the absence of nonlocal multipath, s_(R)(t;ℓ) is modeled by

$\begin{matrix} \begin{array}{l} {s_{R}\left( {t_{0} + t;\mathcal{l}} \right) \approx e^{j2\pi{({\alpha_{TR}{(\mathcal{l})}t + \frac{1}{2}\alpha_{TR}^{(1)}{(\mathcal{l})}t^{2}})}}g_{TR}\left( \mathcal{l} \right)\left( {\text{a}_{TR}\left( {t;\mathcal{l}} \right) \circ s_{T}\left( \left( {1 -} \right) \right)} \right)} \\ {\left( \left( {\left( {\tau_{TR}^{(1)}\left( \mathcal{l} \right)} \right)t - \tau_{TR}\left( \mathcal{l} \right);\mathcal{l}} \right) \right),} \end{array} & \text{­­­Eqn (74)} \end{matrix}$

where a_(TR)(t;ℓ) is the M_(feed) ×1 time-varying spatial signature operator with frequency response A_(TR)(f;ℓ), which is also assumed to adhere to a first-order spatial signature blur model in presence of channel dynamics, due to both movement of the receiver over the reception interval, and adjustments in the yaw, pitch, and roll orientation of the receiver platform. Assuming the antennas have nonzero gain along right-hand and left-hand circular polarizations, and in absence of local scattering multipath and channel dynamics, a_(TR)(t;ℓ) can be modeled as

$\begin{matrix} {\text{A}_{TR}\left( {f;\mathcal{l}} \right) = \text{A}_{R}\left( {f;\Psi_{R}\text{u}_{TR}\left( \mathcal{l} \right)} \right)\rho_{TR}\left( \mathcal{l} \right),} & \text{­­­Eqn (75)} \end{matrix}$

where u_(TR)(ℓ) is the observed line-of-bearing (LOB) direction vector given in Eqn (15) and Ψ_(R) is a 3×3 rotation matrix that captures the yaw, pitch, and tilt of the receiver platform and converts the LOB vector to a local unit-norm direction-of-arrival (DOA) vector, in some aspects parameterized with respect to azimuth and elevation angles, and where

$\left\{ {\text{A}_{R}\left( {f;\text{u}} \right)} \right\}_{{\|\text{u}\|}_{2} = 1} = \left\{ \left\lbrack \begin{array}{ll} \left( {\text{A}_{R}\left( {f;\text{u}} \right)} \right|_{\text{RHCP}} & \left( {\text{A}_{R}\left( {f;\text{u}} \right)} \right|_{LHCP} \end{array} \right\rbrack \right\}_{{\|\text{u}\|}_{2} = 1}$

is the array manifold of M_(feed)×2 complex gains along the right-hand and left-hand circular polarizations at the BFN output, parameterized with respect to 3×1 unit-norm local direction-of-arrival (DOA) vector u, and ρ_(TR)(ℓ) is a 2×1 unit-norm polarization gain vector. The array manifold can include adjustments to account for multipath local to the receiver platform, mutual coupling between antenna elements, and direction-independent complex gains in any distribution system coupling the array to the BFN and/or the receiver bank.

As shown in FIG. 13 , each element of signal vector x_(R)(t) is then channelized into a K_(chn)×1 signal vector

x_(chn)(n_(sym); m_(feed)) = [x_(chn)(k_(chn), n_(sym); m_(feed))]_(k_(chn) = 0)^(K_(chn) − 1),

using a bank of M_(feed) symbol-rate synchronous 1: K_(chn)vector channelizers (250), for example, using M_(feed) 1: M_(ADC) S/P operations and K_(chn)×M_(ADC) channelization matrix operators T_(chn) (253), to form channelized signal x_(chn)(n_(sym);m_(feed)) = T_(chn)◦[x_(R)(n_(sym)T_(sym)+m_(ADC)T_(ADC);m_(feed))]. The Mfeed channelizer output signals

{x_(chn)(n_(sym); m_(feed))}_(m_(feed) = 1)^(M_(feed))

are then combined to form a single M_(chn)×1 complex vector x_(chn)(n_(sym)) (259), where M_(chn)=K_(chn)M_(feed). In some aspects of the invention, this is performed by “stacking” the signals first by receiver feed, and then by channel, such that

$\begin{matrix} {\text{x}_{\text{chn}}\left( n_{\text{sym}} \right) = \begin{pmatrix} \left\lbrack {x_{\text{chn}}\left( {0,n_{\text{sym}};m_{\text{feed}}} \right)} \right\rbrack_{m_{\text{feed}} = 1}^{M_{\text{feed}}} \\  \vdots \\ \left\lbrack {x_{\text{chn}}\left( {K_{\text{chn}} - 1,n_{\text{sym}};m_{\text{feed}}} \right)} \right\rbrack_{m_{\text{feed}} = 1}^{M_{\text{feed}}} \end{pmatrix}.} & \text{­­­Eqn (76)} \end{matrix}$

In others, combining is performed first by channel, then by receiver feed, such that

x_(chn)(n_(sym))=

[x_(chn)(n_(sym); m_(feed))]_(m_(feed) = 1)^(M_(feed)).

Other, more general combining strategies can also be used; however, the form shown in Eqn (76) is has strong advantages for multi-subband processing.

FIG. 14 depicts one aspect for performing multifeed, multi-subband channelization, using an FFT-based frequency channelization method. The M_(feed) ×1 vector ADC output signal x_(ADC)(n_(ADC))= x_(R)(T_(ADC)N_(ADC)) is passed through a bank of M_(feed) 1:M_(ADC) serial-to-parallel (S/P) convertors (251) and optionally windowed (153) symbol-rate synchronous FFT operations (252), each of which computes a K_(FFT) -bin FFT using M_(ADC) contiguous samples of ADC output data, where f_(ADC) = M_(ADC)f_(sym). The FFT operation creates K_(FFT)M_(feed)×1 data vector

[x_(FFT)(k_(bin), n_(sym))]_(k_(bin) = 0)^(K_(FFT) − 1)

with rate f_(sym), where

x_(FFT)(k_(bin); n_(sym)) = [x_(FFT)(k_(bin), n_(sym); m_(feed))]_(m_(feed) = 0)^(M_(feed)).

The FFT output vector is then passed through a subband selection operation (254), which selects K_(sub) subbands of FFT bins to form subbands

{x_(DoF)(n_(sym); k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1),

each with rate f_(sym). In the aspect shown in FIG. 14 , each of the subbands comprise M_(sub) FFT bins, such that x_(DoF)(n_(sym);k_(sub)) is an

$\begin{array}{l} {M_{\text{DoF}} \times 1\text{vector given by x}_{\text{DoF}}\left( {n_{\text{sym}};k_{\text{sub}}} \right) = \left\lbrack {\text{x}_{\text{FFT}}\left( {k_{\text{bin}}\left( {m_{\text{sub}};k_{\text{sub}}} \right),} \right)} \right)} \\ {\left( \left( n_{\text{sym}} \right) \right\rbrack_{m_{\text{sub}} = 0}^{M_{\text{sub}} - 1},} \end{array}$

and where M_(DoF)= M_(sub)M_(feed). Other aspects may vary the number of FFT bins in each subband, e.g., to account for interference loading or spectral shaping issues in each subband. In the aspect shown in FIG. 14 , a “dense” subband selection strategy in which k_(bin)(m_(sub);k_(sub))=(k_(init)(k_(sub))+m_(sub))modK_(FFT) is also shown; however, as in the single-feed receiver, different strategies can be used to form each subband set.

FIG. 15 depicts another aspect of the invention, in which the frequency channelization is performed using a bank of polyphase filter based symbol-rate synchronous channelizers. The M_(feed) ×1 vector ADC output signal x_(ADC)(n_(ADC))= x_(R)(T_(ADC)n_(ADC)) is passed through a bank of M_(feed) 1:M_(ADC) serial-to-parallel (S/P) convertors (251) and a bank of M_(feed) polyphase filter based channelizers (255), each using channelizer weights {w_(chn)(m_(ADC))} (158), which creates a K_(chn)M_(feed)×1 data vector

[x_(chn)(k_(chn), n_(sym))]_(k_(chn) = 0)^(K_(chn) − 1)

with rate f_(sym), where

x_(chn)(k_(chn), n_(sym)) = [x_(chn)(k_(chn), n_(sym); m_(feed))]_(m_(feed) = 0)^(M_(feed)).

As in the single-feed channelizer shown in FIG. 10 , the frequency channels are given by Eqn (21).

The K_(chn)M_(feed)×1 data vector is then put through a M_(DoF)×K_(sub) reshape operation (161), where K_(chn)=K_(sub)M_(sub), and M_(DoF)=M_(sub)M_(feed), such that K_(chn)M_(feed)=K_(sub)M_(DoF). This creates K_(sub) subband signals

{x_(DoF)(n_(sym); k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1),

where M_(DoF)×1 subband k_(sub) signal

x_(DoF)(n_(sym); k_(sub))=

[x_(chn)(K_(sub)k_(sub) + m_(sub), n_(sym))]_(m_(sub) = 0)^(M_(sub) − 1).

FIG. 16 and FIG. 17 show exemplary values of K_(chn)and K_(sub), and M_(sub), designed to yield a constant value of M_(DoF) = 60 as the number of receiver feeds is varied from 1 to 6, for the GPS L1 legacy signal (FIG. 16 ) and the GPS L5 civil signal (FIG. 17 ). In all cases, the receiver channelizer spans roughly 96% of the signal mainlobe, or just past the half-power signal bandwidth. As the number of receiver feeds grows, the number of frequency channels per subband is dropped, and the number of subbands is increased correspondingly, so that the degrees of freedom for linear-algebraic processing is the same in each subband, and the multi-subband channelizer continues to occupy the same overall bandwidth.

The order-of-magnitude complexity in Mop/second and convergence time in milliseconds is also shown in FIG. 16 and FIG. 17 . As these Figures show, the receiver complexity grows linearly with the number of receiver feeds and bandwidth, and is well within the capabilities of low-cost DSP for the L1 legacy signal, and moderate cost DSP for the GPS L5 signal. In contrast, the needed convergence time can be a low multiple of 60 ms for both signals and all receiver feeds. Moreover, the value of FIG. 16 and FIG. 17 chosen here should be sufficient to detect and demodulate all of the GPS signals within the field of view of any terrestrial receiver, with margin left over for detection and removal of pseudolites and spoofers possibly impinging on the GPS band. Lastly, the multifeed receiver should be able excise as many as M_(feed)-1 wideband non-GPS signals, e.g., jammers, or as many as K_(sub) -1 narrowband jammers.

The channel model developed for the single-feed receiver also extends in a straightforward fashion to the multifeed receiver. Ignoring channel dynamics, the polyphase-filter based channelizer yields the same channelized receive signal model given in Eqn (25), and the subband model given in Eqn (58), where

$\begin{matrix} {\text{a}_{\text{chn}}\left( \mathcal{l} \right) \approx \left\lbrack {a_{\text{chn}}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right);\mathcal{l}} \right)\text{A}_{TR}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right);\mathcal{l}} \right)} \right\rbrack,} & \text{­­­Eqn (77)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {\text{a}_{\text{DoF}}\left( {k_{\text{sub}};\mathcal{l}} \right) \approx \left\lbrack {\text{a}_{\text{chn}}\left( {f_{\text{chn}}\left( {K_{\text{sub}} + m_{\text{sub}}} \right);\mathcal{l}} \right)} \right\rbrack_{m_{\text{sub}} = 0}^{M_{\text{sub}} - 1} \otimes} \\ {\text{A}_{TR}\left( {f_{\text{sub}}\left( k_{\text{sub}} \right)\mathcal{l}} \right),} \end{array} & \text{­­­Eqn (78)} \end{matrix}$

and where a_(chn)(f_(chn)(k_(chn));ℓ) is given by Eqn (40)-Eqn (41) and “⊗” denotes the Kronecker product operation. Similarly, the per-channel and per-subband interference ACM is given by

$\begin{matrix} \begin{array}{l} {\text{R}_{\text{i}_{\text{chn}}\text{i}_{\text{chn}}}\left( k_{\text{chn}} \right) \approx \frac{M_{\text{ADC}}^{2}}{M_{\text{chp}}}f_{\text{chp}}\text{diag}} \\ {\left\{ {\left| {H_{R}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)} \right|^{2}\text{S}_{\text{i}_{R}\text{i}_{R}}\left( {f_{\text{chn}}\left( k_{\text{chn}} \right)} \right)} \right\}_{k_{\text{sub}}, = 0}^{K_{\text{sub}} - 1},} \end{array} & \text{­­­Eqn (79)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {\text{R}_{\text{i}_{\text{DoF}}\text{i}_{\text{DoF}}}\left( k_{\text{sub}} \right) \approx \frac{M_{\text{ADC}}^{2}}{M_{\text{chp}}}f_{\text{chp}}\left| {H_{R}\left( {f_{\text{sub}}\left( k_{\text{sub}} \right)} \right)} \right|^{2}} \\ {\left( {\text{P}_{\text{sub}} \otimes \text{S}_{\text{i}_{R}\text{i}_{R}}\left( {f_{\text{sub}}\left( k_{\text{sub}} \right)} \right)} \right),} \end{array} & \text{­­­Eqn (80)} \end{matrix}$

where S_(iRiR) (f) is the M_(feed) ×M_(feed) matrix power spectral density of i_(R)(t).

Importantly, while the TOA and FOA of a GNSS transmitter can be easily spoofed in a covert or “aligned” spoofing scenario, the DOA (and, to a lesser degree, the polarization) of that transmitter cannot be easily spoofed. In addition, the multi-feed receiver can null any CCl impinging on the array, if the array has sufficient degrees of freedom to separate that CCI from the GNSS signals.

In presence of channel dynamics, the local signal DOA adheres closely to a first-order model over time intervals on the order of 10 seconds in some aspects of the invention, and the individual subbands will experience additional signature blur due to the changing TOA, LOB, and (typically more importantly, due to dependence on changing receiver platform orientation) DOA of the transmitters and receiver. This signature blur is likely to further load the subbands with low-level signature components that will be excised by subsequent linear combining operations. In this regard, the effect of TOA changes is reduced substantively for the multifeed symbol-rate synchronous sub-band channelizer shown in FIG. 15 , due to the narrower subband channels obtaining as the number of receiver feeds is increased, as shown in FIG. 16 and FIG. 17 .

RESILIENT PNT SIGNAL PROCESSING METHODS

A common set of signal processing methods can be applied to data provided by all of the single-feed and multifeed receiver structures and symbol-rate synchronous channelization and multi-subband formation operations described above. In particular, in-subband linear combining weights approaching the max-SINR solution can be computed using linear-algebraic methods disclosed in the ‘417 NPA, given any of the following:

-   Knowledge of the content of -   {d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)), -   e.g., from navigation data or symbol sequences provided by     third-party sources. -   Knowledge of a component of -   {d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)), -   e.g., known or periodic pilot signals, as listed in FIG. 4 for the     GPS civil L5, QZSS, Galileo E5AB, E6B E6C, and NAVIC RS BOC signals -   Other exploitable structure of -   {d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)), -   e.g., real or BPSK structure, as listed in FIG. 4 for the GLONASS     and GPS L1 legacy signals; QPSK structure, as listed in FIG. 4 for     the Beidou 1B and NAVIC SPS-L5 and SPS-S signals; and constant     modulus structure, as listed in FIG. 4 for all of the GNSS signals.

These methods also can be used to estimate the max-SINR obtaining in each subband, thereby providing cross-subband weights for the second combining stage shown in Eqn (64). The multi-subband solution can further be used to estimate the full max-SINR of the received symbol sequences, detect the navigation signals in the environment, determine at least the geo-observables

{n_(TR)(𝓁), α̃_(TR)(𝓁)}_(𝓁 = 1)^(L_(T))

for those signals, and if needed estimate information-bearing components of the transmitted symbol streams

{d_(T)(n_(sym); 𝓁)}_(𝓁 = 1)^(L_(T)).

oreover, these techniques can perform these operations without prior knowledge of the signal ranging code or receive frequency signature for the signals, or a search over TOA and FOA is the ranging code is known.

In addition, the ‘417 NPA discloses means for estimating the fine TOA, FOA Nyquist zone, and DFOA using additional copy-aided parameter estimation algorithms. Lastly, the data model derived above can be used to develop matched-filtering methods, using estimates of a_(chn)(ℓ) given in Eqn (40)-Eqn (41) for single-receiver systems, or defined in Eqn (77) for multi-receiver systems.

The signal processing structures can be adapted on either a continuous basis, in which geo-observables are updated rapidly over time, or on a batch processing basis, in which a block of N_(sym) channelized data symbols

{x_(DoF)(n_(sym); k_(sub))}_(n_(sym) = 0)^(N_(sym) − 1)

are computed for each subband and passed to a DSP processing element that detects the GNSS signals within that data block. The latter approach is especially useful if the invention is being used to develop resilient PNT analytics to aid a primary navigation system, e.g., to assess quality and availability of new GNSS transmissions, or to detect or confirm spoofing transmissions on a periodic basis. The batch adaptive processing algorithms are described in more detail below.

Batch Adaptive Processing Procedure

In the batch adaptive processing procedure is implemented by first collecting data over N_(sym) data symbols, and detecting, extracting, or estimating geo-observables of the MOS-DSSS symbol or navigation sequences directly from that data set. In some aspects, the procedure performs this processing from a “cold start,” i.e., with no prior information about the signals contained within that data set. However, if the prior FOA’s of the signals (and in particular FOA’s derived from the FOA vectors for those signals) are known, then the procedure can be started at an intermediate point in the processing.

This procedure enables a great deal of refinement and accurate discrimination to more closely constrain and limit the processing necessary to accurately interpret the signal’s content, before the copy-aided analysis phase begins. In some use scenarios, the blind despreading stage can in fact obviate the copy-aided analysis phase, e.g., if the invention is developing resilient PNT analytics to aid a primary navigation system, or it can be used to substantively thin the number of transmissions that must be analyzed. This procedure can thereby reduce the processing complexity and considerable feedback lag, enabling quicker, more effective signal discrimination without requiring the full processing and analysis of the signal be completed first (or even together).

FIG. 18 is a flow diagram of a fully-blind multi-subband signal detection, geo-observable/quality estimation, and symbol estimation batch procedure that can be employed for MOS-DSSS signals with BPSK (or, more generally, real) symbol sequences, e.g., the GLONASS and GPS L1 legacy signal, or LocataLite beacons. Upon receipt of N_(sym) symbols of data (300), e.g., N_(sym) milliseconds of receive data for the GNSS signals depicted in FIG. 4 , the procedure first channelizes the data and forms it into K_(sub) subbands using a symbol-rate synchronous channelizer (301). In some aspects, the data subbands are formed into K_(sub) data matrices with dimension N_(sym) ×M_(DoF), e.g., given by

$\begin{matrix} {\text{x}_{\text{DoF}}\left( k_{\text{sub}} \right) \triangleq \begin{pmatrix} {\sqrt{\omega_{\text{DoF}}(0)}\text{x}_{\text{DoF}}^{H}\left( {0;k_{\text{sub}}} \right)} \\  \vdots \\ {\sqrt{\omega_{\text{DoF}}\left( {N_{\text{sym}} - 1} \right)}\text{x}_{\text{DoF}}^{H}\left( {N_{\text{sym}} - 1;k_{\text{sub}}} \right)} \end{pmatrix},} & \text{­­­Eqn (81)} \end{matrix}$

where

{ω_(DoF)(n_(sym))}_(n_(sym) = 0)^(N_(sym) − 1)

is a real data window satisfying

$\sum\limits_{n_{\text{sym}} = 0}^{N_{\text{sym}} - 1}{\omega_{\text{DoF}}\left( n_{\text{sym}} \right) = 1.}$

The channelized subband data matrices

{X_(DoF)(k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1)

are then whitened (302), e.g., by performing the QR decomposition (QRD) of each subband matrix, to form whitened subband matrices

{Q_(DoF)(k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1).

The QRD, denoted {Q,R} = QRD(X) for general N×M matrix X, solves

$\begin{matrix} {\text{R}\text{=}\mspace{6mu}\mspace{6mu}\text{chol}\left\{ {\text{X}^{H}\text{X}} \right\},} & \text{­­­Eqn (82)} \end{matrix}$

$\begin{matrix} {\text{Q}\mspace{6mu}\mspace{6mu}\text{=}\mspace{6mu}\text{XR}^{- 1} = \begin{pmatrix} {\text{q}^{H}(0)} \\  \vdots \\ {\text{q}^{H}\left( {N - 1} \right)} \end{pmatrix},} & \text{­­­Eqn (83)} \end{matrix}$

where chol(↺) is the Cholesky factorization operation, such that X = QR and Q^(H)Q = I_(M)

where I_(M)

M×M identity matrix. The data window is nominally rectangular; however, other windows are recommended if the FOA offset between the navigation signals is small, e.g., for fixed ground beacons or aligned spoofers, or for reception intervals that are long enough to induce substantive differential FOA effects.

For navigation signals in which only a single rail is modulated, e.g., GLONASS or GPS L1 legacy signals, the subbands can be processed using the conjugate self-coherence restoral (CSCORE) method disclosed in the ‘417 NPA. For the multi-subband method, CSCORE statistics are developed within each subband, and combined to form a cross-subband statistic (303) that simultaneously detects the signals, and provides an estimate the fine FOA and quality of each signal. In one aspect, the CSCORE algorithm is implemented in each subband for a trial FOA vector α= [ α α⁽¹⁾]^(T), by first computing

Q_(DoF)(α;k_(sub))  =  Δ(α)Q_(DoF)(k_(sub)),

Δ(α)  ≜  diag{exp(j2πg_(FOA)(n_(sym))α)}_(n_(sym) = 0)^(N_(sym) − 1),

where

$\text{g}_{\text{FOA}}\left( n_{\text{sym}} \right) = \left\lbrack {n_{\text{sym}}\mspace{6mu}\mspace{6mu}\frac{1}{2}n_{\text{sym}}^{2}} \right\rbrack;$

forming whitened CSCORE matrix

${\overset{˙}{R}}_{q_{DoF}q_{DoF}^{*}}\left( {\alpha;k_{sub}} \right)$

given by

$\begin{matrix} \begin{matrix} {{\hat{\text{R}}}_{\text{q}_{\text{DoF}}\text{q}_{\text{DoF}}^{\ast}}\left( {\text{α;}k_{\text{sub}}} \right) = {\sum\limits_{n_{\text{sym}} = 0}^{N_{\text{sym}} - 1}{\text{q}_{\text{DoF}}\left( {n_{\text{sym}};k_{\text{sym}}} \right)\text{q}_{\text{DoF}}^{T}\left( {n_{\text{sym}};k_{\text{sym}}} \right)}}} \\ {\exp\left( {- j4\pi\text{g}_{\text{FOA}}^{T}\left( n_{\text{sym}} \right)\text{α}} \right)} \\ {= {\sum\limits_{n_{\text{sym}} = 0}^{N_{\text{sym}} - 1}\text{q}_{\text{DoF}}}\left( {n_{\text{sym}};\text{α;}k_{\text{sub}}} \right)\text{q}_{\text{DoF}}^{T}\left( {n_{\text{sym}};\text{α;}k_{\text{sub}}} \right),} \\ {= \text{Q}_{\text{DoF}}^{H}\left( {\text{α;}k_{\text{sub}}} \right)\text{Q}_{\text{DoF}}^{\ast}\left( {\text{α;}k_{\text{sub}}} \right)} \end{matrix} & \text{­­­Eqn (84)} \end{matrix}$

for each subband; and determining the dominant mode {λ_(CSC)(α;k_(sub)),v_(CSC)(α;k_(sub))} of the singular value decomposition (SVD) of

${\overset{˙}{R}}_{q_{DoF}q_{DoF}^{*}}\left( {\alpha;k_{sub}} \right),$

e.g., using a power method. The detection statistics

{λ_(CSC)(α;k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1)

are then combined across the subbands to form a full-band CSCORE statistic. In one aspect on the invention, this is accomplished by first computing max-SINR estimate

γ_(CSC)(α;k_(sub))=

(1 + λ_(CSC)(α;k_(sub)))/(1 − λ_(CSC)(α;k_(sub)))

for each subband, and summing those estimates together foreach FOA vector α. In other aspects, nonlinear combining operations, e.g., dictated by maximum-likelihood (ML) estimation arguments are performed. In other aspects, the combining is performed over values of α that are adjusted to account for DTOA between subbands.

This procedure generates a CSCORE spectrum that is parameterized with respect to both FOA (or DTOA) and DFOA. This spectrum is then used to detect the MOS-DSSS signals in the environment, and estimate their FOA (or DTOA) and DFOA geo-observables (304). Additional joint processing is then performed to further improve quality of the geo-observables and the detection statistic (305), resulting in optimized DTOA and DFOA estimates

{τ̂_(CSC)⁽¹⁾(𝓁), α̂_(CSC)⁽¹⁾(𝓁)}_(𝓁 = 1)^(L̂_(T))

for the L _(T) detected signals, and resulting in optimized despreading weights

{v_(CSC)(α(k_(sub); 𝓁)); k_(sub)}_(𝓁 = 1)^(L̂_(T))

for each subband, where

α(k_(sub); 𝓁)=

[−τ̂_(CSC)⁽¹⁾(𝓁)(f_(T) + f_(sub)(k_(sub)))  α̂_(CSC)⁽¹⁾(𝓁)]^(T).

The optimized despreading weights and estimated FOA vector are then used to despread and despin the symbol sequence for each detected signal, and to demodulate the underlying signal navigation sequence based on further structure of that sequence, e.g., 20:1 NAV bit replication for the GPS L1 legacy signal (306). As part of this process, the coarse TOA timing is computed to within a NAV bit edge.

Given reception of sufficient data, the NAV signal is analyzed to resolve the ±1 ambiguity in the spread signal; detect NAV block boundaries to determine the full coarse TOA; and if needed extract the satellite ephemeres from the NAV signal sequences (307). In some aspects, this is accomplished by processing of data over multiple reception blocks. In other aspects, this is accomplished by shifting to a continuous update (non-batch) processing mode.

In some aspects, the ranging code and (for multifeed receivers) array manifold are downloaded from memory (310), and used to determine the full TOA, FOA, DOA, and (using estimated and/or on-board orientation data) LOB of the detected signals, and determine a positioning, navigation, and timing (PNT) solution for the receiver (308). In other aspects, the fine FOA (or TDOA) and DFOA, and transmitter ephemeres are sufficient to determine a positioning, navigation and frequency synchronization solution for the receiver.

FIG. 19 is a flow diagram of a multi-subband signal detection and geo-observable/quality estimation, algorithm that can be employed in a resilient PNT data analytics engine (RPNT-AE) aspect of the invention. Upon receipt of an RPNT-AE trigger (319), e.g., initiated by a primary GNSS processor, or by other processes on board and/or connected to the receiver, the RPNT-AE process is begun (320). The RPNT-AE processor then tasks a receiver to obtain and store a snapshot of N_(sym) data symbols, and obtains a time/frequency stamp recording rough (to within receiver clock accuracy) center frequency and start time of the snapshot, and other information available to the receiver, e.g., receiver orientation (321). The RPNT-AE processor then places the time/frequency stamped data into memory (310), and sends the stamp to an external resource (322 a) along with a request for third-party data. The RPNT-AE processor then waits until an external resource (322 b) provides the processor with a snapshot of navigation (NAV/CNAV) data covering the time period and frequency band recorded in the time/frequency stamp, and (if needed) a transmitter almanac containing the ephemeres of the transmitters over the time period covering the NAV/CNAV data. Once RPNT-AE processor obtains that data (323), it places that data in memory (310) and sends a “Data Ready” flag is sent to the DSP processor tasked with estimating RPNT analytics for the processor.

Once the Data Ready flag is obtained from the RPNT-AE processor, the DSP resources obtain the time/frequency stamped received data snapshot and the NAV/CNAV data from memory (310), and generate RPNT analytics using partially blind methods disclosed in the ‘417 NPA, e.g., an “FFT-least-squares” algorithm or a single-target or multi-target maximum-likelihood estimator, the DSP processor. Analytics measure here can include the following:

-   estimates of observed received parameters of navigation signals,     including signal geo-observables usable to obtain PNT solutions for     the reception platform, e.g., observed signal frequency-of-arrival     (FOA), time-of-arrival (TOA), and observed local     direction-of-arrival (DOA) or global line-of-bearing (LOB) in     systems employing multifeed receivers; -   estimates of received signal quality, e.g., received incident power,     and of despread/demodulated navigation sequence quality, e.g.,     despreader output signal-to-interference-and-noise ratio (SINR); and -   measures of accuracy of parameter and signal quality estimates.

RPNT analytic measurement methods include “partially-blind” methods that do not require knowledge of the ranging code for the navigation signals, or spatial/polarization array-manifold data (measurement of cross-feed spatial/polarization signatures as a function of DOA) for the receiver, and “copy-aided” methods that may require one or both of the ranging codes or the array manifold to provide more complete RPNT analytics. In the aspects disclosed herein, RPNT analytics are further computed over multiple subbands.

Once the RPNT analytics have been obtained, the RPNT-AE process is ended (325), and the RPNT analytics are reported to the resources requesting the analytics.

FIG. 20 is a block diagram of a multi-subband pilot-exploiting RPNT (PE-RPNT) processor that can detect and estimate geo-observables and quality of MOS-DSSS signals in which known or repeating pilots transmitted as part of the MOS-DSSS system sequence. The method reported herein can also be used to implement the operations performed by the RPNT-AE DSP processor (324), using the full MOS-DSSS symbol sequence provided in that system.

Upon receipt of N_(sym) symbols of data (300), e.g., N_(sym) milliseconds of receive data for the GNSS signals depicted in FIG. 4 , the procedure first channelizes the data and forms it into K_(sub) subbands using a symbol-rate synchronous channelizer (301). In some aspects, the data subbands are formed into K_(sub) data matrices with dimension N_(sym) ×M_(DoF), e.g., given by Eqn (81). The channelized subband data matrices

{X_(DoF)(k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1)

are then whitened (302), e.g., by performing the QR decomposition (QRD) of each subband matrix, to form whitened subband matrices

{Q_(DoF)(k_(sub))}_(k_(sub) = 0)^(K_(sub) − 1).

The processor then retrieves a single period of repeating pilot data

{p(n_(sym))}_(n_(sym) = 0)^(N_(pilot) − 1)

from memory (310), where N_(pilot) is the period of the pilot sequence; computes pilot statistics over a set of trial FOA vectors α=[ α α⁽¹⁾ ]^(T) trial pilot delay m_(pilot) ∈ {0,,N_(pilot)-1} in each subband, given by least-squares (LS) whitened linear combiner weights

$\begin{matrix} \begin{matrix} {\text{u}_{\text{LS}}\left( {\text{α,}m_{\text{pilot}};k_{\text{sub}}} \right) = {\sum\limits_{n_{\text{sym}} = 0}^{N_{\text{sym}}}{\text{q}_{\text{DoF}}\left( {n_{\text{sym}};k_{\text{sub}}} \right)d_{\text{LS}}^{\ast}\left( {n_{\text{sym}};\text{α,}m_{\text{pilot}}} \right)}}} \\ {= \text{Q}_{\text{DoF}}^{H}\left( k_{\text{sub}} \right)\text{d}_{\text{LS}}\left( {\text{α,}m_{\text{pilot}}} \right),} \end{matrix} & \text{­­­Eqn (85)} \end{matrix}$

$\begin{matrix} \begin{matrix} {d_{\text{LS}}\left( {n_{\text{sym}};\text{α,}m_{\text{sym}}} \right) = \sqrt{\omega_{\text{DoF}}\left( n_{\text{sym}} \right)}p\left( {\left( {n_{\text{sym}} - m_{\text{pilot}}} \right)\text{mod}N_{\text{pilot}}} \right)} \\ e^{j2\pi\text{g}_{\text{FOA}}{(n_{\text{sym}})}\text{α}} \\ {\text{d}_{\text{LS}}\left( {\text{α,}m_{\text{sym}}} \right) = \begin{pmatrix} {d_{\text{LS}}^{\ast}\left( {0;\text{α,}m_{\text{sym}}} \right)} \\  \vdots \\ {d_{\text{LS}}^{\ast}\left( {N_{\text{sym}} - 1;\text{α,}m_{\text{sym}}} \right)} \end{pmatrix},} \end{matrix} & \text{­­­Eqn (86)} \end{matrix}$

where

$\text{g}_{\text{FOA}}\left( n_{\text{sym}} \right) = \left\lbrack {n_{\text{sym}}\mspace{6mu}\mspace{6mu}\frac{1}{2}n_{\text{sym}}^{2}} \right\rbrack,$

and in-band LS SINR estimate,

$\begin{matrix} {\text{γ}_{\text{LS}}\left( {\text{α},m_{\text{pilot}};k_{\text{sub}}} \right) = \frac{\left\| {\text{u}_{\text{LS}}\left( {\text{α},m_{\text{pilot}};k_{\text{sub}}} \right)} \right\|_{2}^{2}}{1 - \left\| {\text{u}_{\text{LS}}\left( {\text{α},m_{\text{pilot}};k_{\text{sub}}} \right)} \right\|_{2}^{2}};} & \text{­­­Eqn (87)} \end{matrix}$

and combines the in-subband least-square SINR estimates to created a multi-subband quality statistic as a function of FOA vector and pilot delay (333). In one aspect, the in-subband LS SINR estimates are combined in accordance with a maximum-likelihood estimator, yielding multi-subband quality statistic

$\begin{matrix} {S_{\text{ML}}\left( {\text{α},m_{\text{pilot}}} \right) = {\sum\limits_{k_{\text{sub}} = 0}^{K_{\text{sub}} - 1}{\log\left( {1 + \text{γ}_{\text{LS}}\left( {\text{α},m_{\text{pilot}};k_{\text{sub}}} \right)} \right)}},} & \text{­­­Eqn (88)} \end{matrix}$

which is maximized at modulo- N_(pilot) coarse TOA’s and fine FOA vectors possessed by the MOS-DSSS signals in the receiver’s field of view. In another aspect, the in-subband LS SINR estimates are combined in accordance with Eqn (65), yielding multi-subband ML quality statistic

$\begin{matrix} {S_{\text{LS}}\left( {\alpha,m_{\text{pilot}}} \right) = {\sum\limits_{k_{\text{sub}} = 0}^{K_{\text{sub}} - 1}{\gamma_{\text{LS}}\left( {\alpha,m_{\text{pilot}};k_{\text{sub}}} \right)}},} & \text{­­­Eqn (89)} \end{matrix}$

which is also maximized at modulo- N_(pilot) coarse TOA’s and fine FOA vectors possessed by the MOS-DSSS signals in the receiver’s field of view. In other aspects, the effect of DTOA across frequency subbands is taken into account, e.g., resulting in DTOA-corrected multi-subband LS quality statistic

$\begin{matrix} {S_{\text{LS}}\left( {\text{τ}^{({1,2})},m_{\text{pilot}}} \right) = {\sum\limits_{k_{\text{sub}}}^{K_{\text{sub}} - 1}{\text{γ}_{\text{LS}}\left( {\text{α}\left( k_{\text{sub}} \right),m_{\text{pilot}};k_{\text{sub}}} \right)}}} & \text{­­­Eqn (90)} \end{matrix}$

$\begin{matrix} {\text{τ}^{({1,2})} = \begin{pmatrix} \text{τ}^{(1)} \\ {\widetilde{\text{τ}}}^{(2)} \end{pmatrix},} & \text{­­­Eqn (91)} \end{matrix}$

$\begin{matrix} {\text{α}\left( k_{\text{sub}} \right) = \mspace{6mu} - \left( {f_{T} + f_{\text{sub}}\left( k_{\text{sub}} \right)} \right)T_{\text{sym}}\text{τ}^{({1,2})},} & \text{­­­Eqn (92)} \end{matrix}$

where τ̃⁽²⁾= τ⁽²⁾T_(sym) is the symbol-normalized differential DTOA. In addition, for rectangular time windows, the in-band LS SINR estimate can be converted to unbiased max-SINR estimate

$\begin{matrix} \begin{array}{l} {{\hat{\gamma}}_{\text{max-SINR}}\mspace{6mu}\left( {\text{α},m_{\text{pilot}};k_{\text{sub}}} \right) = \left( {1 - \frac{M_{\text{DoF}} + 1}{N_{\text{sym}}}} \right)\gamma_{\text{LS}}\left( {\text{α},m_{\text{pilot}};k_{\text{sub}}} \right) -} \\ {M_{\text{DoF}},} \end{array} & \text{­­­Eqn (93)} \end{matrix}$

which can greatly improve visible quality of S_(LS)(α,m_(pilot)) (or S_(LS)(τ^((1,2));m_(pilot))) without affecting performance of the statistic.

A number of means can be used to minimize efficiency of the in-subband LS SINR’s and multi-subband quality statistics. In particular, FFT-based methods can be used to efficiently compute γ_(LS)(α,m_(pilot);k_(sub)) or γ_(LS)(α(k_(sub)),m_(pilot);k_(sub)) with fine accuracy. At this stage of processing, the whitened LS weights given in Eqn (85) need not be computed, further reducing complexity and memory requirements of the overall processor.

The multi-subband quality statistic is then analyzed to detect the pilot-bearing MOS-DSSS signals in the environment, determine their coarse TOA to modulo- N_(pilot) accuracy, and determine their fine FOA or DTOA vector (334), resulting in L̂_(T) peak detections and FOA-TOA locations

{α(𝓁), m_(pilot)(𝓁)}_(𝓁=1)^(L̂_(T)) or DTOA-

TOA locations

{τ^((1, 2))(𝓁), m_(pilot)(𝓁)}_(𝓁=1)^(L̂_(T))

of those peaks. Iterative refinement of the DTOA-TOA locations

(where the FOA vector is converted to a DTOA vector) is then performed (335), e.g., using Newton search methods disclosed in the ‘417 NPA, to determine the fine DTOA vector (from which FOA can be derived) and optimize the multi-subband quality statistic.

Using the optimized DTOA vectors and modulo-N_(pilot) coarse TOA’s

{τ^((1, 2))(𝓁), m_(pilot)(𝓁)}_(𝓁=1)^(L_(T)), LS

combiner weights {γ_(LS)(k_(sub);ℓ),u_(LS)(k_(sub);ℓ)} are then computed using formula

$\begin{matrix} \begin{array}{l} {\text{u}_{\text{LS}}\left( {k_{\text{sub}};\mathcal{l}} \right) =} \\ {{\sum\limits_{n_{\text{sym}}\text{=0}}^{N_{\text{sym}}}{\text{q}_{\text{DoF}}\left( {n_{\text{sym}};k_{\text{sub}}} \right)d_{\text{LS}}^{*}\left( {n_{\text{sym}};\alpha\left( {k_{\text{sub}};\mathcal{l}} \right),m_{\text{Pilot}}\left( \mathcal{l} \right)} \right)}},} \end{array} & \text{­­­Eqn (94)} \end{matrix}$

$\begin{matrix} {\text{α}\left( {k_{\text{sub}};\mathcal{l}} \right) = - \left( {f_{T} + f_{\text{sub}}\left( k_{\text{sub}} \right)} \right)T_{\text{sym}}\text{τ}^{({1,2})}\left( \mathcal{l} \right),} & \text{­­­Eqn (95)} \end{matrix}$

$\begin{matrix} {\gamma_{\text{LS}}\left( {k_{\text{sub}};\mathcal{l}} \right) = \frac{\left\| {\text{u}_{\text{LS}}\left( {k_{\text{sub}};\mathcal{l}} \right)} \right\|_{2}^{2}}{1 - \left\| {\text{u}_{\text{LS}}\left( {k_{\text{sub}};\mathcal{l}} \right)} \right\|_{2}^{2}},} & \text{­­­Eqn (96)} \end{matrix}$

and cross-subband combiner gains are computed based on Eqn (64), i.e., by setting g_(LS)(k_(sub);ℓ) = 1+γ_(LS)(k_(sub);ℓ). The full multi-subband despread, despun and partially time-synchronized symbol sequence is then estimated for each peak using formula

$\begin{matrix} {{\hat{d}}_{T}\left( {n_{\text{sym}};\mathcal{l}} \right) = {\sum\limits_{k_{\text{sub}}\text{=0}}^{K_{\text{sub}} - 1}{g_{\text{LS}}^{*}\left( {k_{\text{sub}};\mathcal{l}} \right){\hat{d}}_{T}\left( {n_{\text{sym}};k_{\text{sub}};\mathcal{l}} \right)}},} & \text{­­­Eqn (97)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {{\hat{d}}_{\text{T}}\left( {n_{\text{sym}};k_{\text{sub}};\mathcal{l}} \right) =} \\ {\left( {\text{u}_{\text{LS}}^{H}\left( k_{\text{sub}} \right)q_{\text{DoF}}\left( {n_{\text{sym}} + m_{\text{pilot}}\left( \mathcal{l} \right);} \right)} \right)e^{- j2\pi\text{g}_{\text{FOA}}{(n_{\text{sym}})}\alpha{({k_{\text{sub}};\mathcal{l}})}}.} \end{array} & \text{­­­Eqn (98)} \end{matrix}$

The navigation-bearing component of the symbol sequence (CNAV for the GPS L5 signal) is then demodulated using information about the structure of the navigation data (346).

Given reception of sufficient data, the CNAV signal is analyzed to determine the full coarse TOA, and if needed extract the satellite ephemeres from the NAV signal sequences (347). In some aspects, this is accomplished by processing of data over multiple reception blocks. In other aspects, this is accomplished by shifting to a continuous update (non-batch) processing mode.

In some aspects, the ranging code and (for multifeed receivers) array manifold are downloaded from memory (310), and used to determine the full TOA, FOA, DOA, and (using estimated and/or on-board orientation data) LOB of the detected signals, and determine a positioning, navigation, and timing (PNT) solution for the receiver (308). In other aspects, the fine FOA (or TDOA) and DFOA, and transmitter ephemeres are sufficient to determine a positioning, navigation and frequency synchronization solution for the receiver.

FIG. 21 is a flow diagram of a local-maxima search procedure according to an aspect of the disclosure, which is applicable to searches over FOA and DFOA, over reception intervals in which DFOA can have a substantive effect on detection sensitivity. The approach assumes the use of an FFT-based coarse search procedure, in which α is defined over search grid

$\begin{matrix} {\widetilde{\alpha}\left( k_{\text{bin}} \right) = {k_{\text{bin}}/K_{\text{bin}}},\quad k_{\text{bin}} = 0,\ldots,K_{\text{bin}} - 1,} & \text{­­­Eqn (99)} \end{matrix}$

$\begin{matrix} {{\widetilde{\alpha}}_{TR}^{(1)}\left( k_{\text{tile}} \right) = \left( {\widetilde{\alpha}}_{TR}^{(1)} \right|_{\max}\left( \frac{2k_{\text{tile}} + 1}{K_{\text{tile}}} \right),\quad k_{\text{tile}} = 0,\ldots,K_{\text{tile}} - 1,} & \text{­­­Eqn (100)} \end{matrix}$

where K_(bin) is the number of FOA bins searched over, and K_(tile) is the number of DFOA tiles searched over. Then the multi-subband quality estimate can be given by 3-dimensional surface

$\begin{matrix} {S_{\text{LS}}\left( {k_{\text{bin}},k_{\text{tile}},m_{\text{pilot}}} \right) = {\sum\limits_{k_{\text{sub}} = 0}^{K_{\text{sub}} - 1}{\gamma_{\text{LS}}\left( {k_{\text{bin}},k_{\text{tile}},m_{\text{pilot}};k_{\text{sub}}} \right)}},} & \text{­­­Eqn (101)} \end{matrix}$

$\begin{matrix} {\gamma_{\text{LS}}\left( {k_{\text{bin}},k_{\text{tile}},m_{\text{pilot}};k_{\text{sub}}} \right) = \frac{\left\| {\text{u}_{\text{LS}}\left( {k_{\text{bin}},k_{\text{tile}},m_{\text{pilot}};k_{\text{sub}}} \right)} \right\|_{2}^{2}}{1 - \left\| {\text{u}_{\text{LS}}\left( {k_{\text{bin}},k_{\text{tile}},m_{\text{pilot}};k_{\text{sub}}} \right)} \right\|_{2}^{2}},} & \text{­­­Eqn (102)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {\text{u}_{\text{LS}}\left( {k_{\text{bin}},k_{\text{tile}},m_{\text{pilot}};k_{\text{sub}}} \right) =} \\ {{\sum\limits_{n_{\text{sym}} = 0}^{N_{sym} - 1}{\text{q}_{\text{DoF}}\left( {n_{\text{sym}};k_{\text{sub}}} \right)d_{\text{LS}}^{*}\left( {n_{\text{sym}};k_{\text{tile}},m_{\text{pilot}}} \right)e^{{- j2\pi k_{\text{bin}}n_{\text{sym}}}/K_{\text{bin}}}}} =} \\ {\text{DFT}\left\{ {\text{q}_{\text{DoF}}\left( {n_{\text{sym}};k_{\text{sub}}} \right)d_{\text{LS}}^{*}\left( {n_{sym};k_{\text{tile}},m_{\text{pilot}}} \right)} \right\},} \end{array} & \text{­­­Eqn (103)} \end{matrix}$

$\begin{matrix} \begin{array}{l} {d_{\text{LS}}\left( {n_{\text{sym}};k_{\text{tile}},m_{\text{pilot}}} \right) =} \\ {e^{j\pi{\widetilde{\alpha}}^{(1)}{(k_{\text{tile}})}n_{\text{sym}}^{2}}p\left( {\left( {n_{\text{sym}} - m_{\text{pilot}}} \right){mod}N_{\text{pilot}}} \right),} \end{array} & \text{­­­Eqn (104)} \end{matrix}$

which can be implemented using K_(tile)N_(pilot) DFT operations.

Returning to FIG. 21 , once the search is started (350), the number of DFOA search tiles K_(tile) first is determined based on the expected range in DFOA and the reception time. Then an initial DFOA search tile is chosen (351), and a DFT-based search is conducted over that tile using Eqn (101)-Eqn (104), to compute a multi-subband search spectrum for that tile, and determine peaks in that tile (352).

Once all of the DFOA tiles have been searched, and peaks have been determined, the detected peaks are thinned to a significant number of peaks based on peak value (353). Then the peak values are iterivaley optimized using local optimization methods, e.g., quadratic fit to the peak maximum, followed by Newton or Gauss-Newton search methods (354). In this case, m_(pilot) is defined over an integer set of values, hence it is not iteratively optimized using this procedure.

In one aspect of the invention, the local optimization is performed to optimize

$\begin{matrix} {S_{\text{LS}}\left( \mathcal{l} \right) = {\sum\limits_{k_{\text{sub}} = 0}^{K_{\text{sub}} - 1}{\gamma_{\text{LS}}\left( {k_{\text{sub}};\mathcal{l}} \right)}},} & \text{­­­Eqn (105)} \end{matrix}$

where γ_(LS)(k_(sub);ℓ) is given by Eqn (94)-Eqn (96), and where

{τ^((1, 2))(𝓁), m_(pilot)(𝓁)}_(𝓁 = 1)^(L̂_(T))

is the parameter set being optimized. In this case, τ^((1,2))(ℓ) is initialized using the optimal value of α(ℓ) determined during the FFT-based search procedure.

Once optimal values have been determined, ancillary RPNT statistics are computed at those optimal search locations (355). The search is then ended (356).

FIGS. 22A, 22B, and 22C, and FIGS. 23A, 23B, and 23C depict DFOA search sensitivity as a function of the number of DFOA search tiles, for a ±3 Hz/s DFOA range, consistent with DFOA’s expected for a MEO satellite constellation. FIGS. 22A, 22B, and 22C show search sensitivity relative to 0 Hz/s for a rectangular time window, and for a ½ second (FIG. 22A), 1 second (FIG. 22B), and 2 second (FIG. 22C) reception interval. As shown in FIGS. 22A, 22B, and 22C, the number of search tiles needed to minimize sensitivity loss grows rapidly with reception time. FIGS. 23A, 23B, and 23C show equivalent search sensitivity for a Nuttall time window, which has a time-bandwidth product (TBP) that is roughly half the value of the rectangular window. The Nuttall window can reduce the number of tiles required to perform an effective search at a 2 second reception time. However, this is offset by the lower TBP of the window.

While this invention is susceptible of instantiation in many different forms, there are shown in the drawings and described in detail in the text of Provisional Appl. No. 62/773,589, incorporated herein by reference, and Provisional Appl. No. 62/773,605, incorporated herein by reference, several specific aspects of the invention, with the understanding that the present disclosure is to be considered as an exemplification of the principles of the invention and is not intended to limit the invention to the aspects illustrated.

In an aspects applicable to all of the approaches above, the invention obtains snapshots of baseband navigation data covering the time interval of data collected by the receiver, and symbol-synchronously channelized by the invention, and uses that baseband navigation data to implement non-fully-blind demodulation algorithms. The resultant aspects can provide extreme high precision of FOA, TOA, and DFOA/DTOA drift estimates; assess integrity of data collected by the host platform; and provide other functions of interest to the user. When coupled with a communication channel allowing the receive data to be transported to a central processor, the invention can also allow implementation of all functions at off-line resources, thereby eliminating all DSP complexity associated with the algorithms. The aspects can also be used to implement signal cancellation algorithms that detect signals under the known navigation signals, e.g., for purposes of spoofer and jammer detection.

Any reception operation used in the invention can be implemented using any of the set of one or more dedicated receivers and software defined radios (SDR) either separate from or integrated with antennas, amplifiers, mixers, filters, analog-to-digital converters (ADC’s) and signal processing gear.

Operations Processing used in each of the inventions above can be implemented in any combination of hardware and software, from special-purpose hardware including any of application -specific integrated circuits (ASIC’s) and field-programmable gate arrays (FPGA’s); firmware instructions in a lesser-specialized set of hardware; embedded digital signal processors (e.g. Texas Instrument or Advanced Risc Machine DSP’s); graphical processing units (GPU’s); vector, polynomic, quantum, and other processors; and in any combination or sole use of serial or parallel processing; and on general-purpose computers using software instructions.

Operations Processing used in each of the inventions above can be further implemented using any set of resources that are on-board, locally accessible to, and remotely accessible by the receiver after transport of the data and instructions to be processed by any of a single computer, server, and set of servers, and then directed onwards, using any number of wired or wireless means for such transport.

Some of the above-described functions may be composed of instructions, or depend upon and use data, that are stored on storage media (e.g., computer-readable medium). Some of the above-described functions may be comprised in EEPROMs, ASICs, or other combinations of digital circuitry for digital signal processing, connecting and operating with the adaptive processor. The instructions and/or data may be retrieved and executed by the adaptive processor. Some examples of storage media are memory devices, tapes, disks, and the like. The instructions are operational when executed by the adaptive processor to direct the adaptive processor to operate in accord with the invention; and the data is used when it forms part of any instruction or result therefrom.

The terms “computer-readable storage medium” and “computer-readable storage media” as used herein refer to any medium or media that participate in providing instructions to a CPU for execution. Such media can take many forms, including, but not limited to, non-volatile (also known as ‘static’ or ‘long-term’) media, volatile media and transmission media. Non-volatile media include, for example, one or more optical or magnetic disks, such as a fixed disk, or a hard drive. Volatile media include dynamic memory, such as system RAM or transmission or bus ‘buffers’. Common forms of computer-readable media include, for example, a floppy disk, a flexible disk, a hard disk, magnetic tape, any other magnetic medium, a CD-ROM disk, digital video disk (DVD), any other optical medium, any other physical medium with patterns of marks or holes.

Memory, as used herein when referencing to computers, is the functional hardware that for the period of use retains a specific structure which can be and is used by the computer to represent the coding, whether data or instruction, which the computer uses to perform its function. Memory thus can be volatile or static, and be any of a RAM, a PROM, an EPROM, an EEPROM, a FLASH-EPROM, any other memory chip or cartridge, a carrier wave, or any other medium from which a computer can read data, instructions, or both.

I/O, or ‘input/output’, is any means whereby the computer can exchange information with the world external to the computer. This can include a wired, wireless, acoustic, infrared, or other communications link (including specifically voice or data telephony); a keyboard, tablet, camera, video input, audio input, pen, or other sensor; and a display (2D or 3D, plasma, LED, CRT, tactile, or audio). That which allows another device, or a human, to interact with and exchange data with, or control and command, a computer, is an I/O device, without which any computer (or human) is essentially in a solipsistic state.

While this invention has been described in reference to illustrative aspects of the invenion, this description is not to be construed in a limiting sense. Various modifications and combinations of the illustrative aspects as well as other aspects of the invention will be apparent to those skilled in the art upon referencing this disclosure. It is therefore intended this disclosure encompass any such modifications or aspects.

The scope of this invention includes any combination of the elements from different aspects disclosed in this specification, and is not limited to the specifics of any of the aspects mentioned above. Individual user configurations and aspects of this invention may contain all, or less than all, of the elements disclosed in the specification according to the needs and desires of that user. The claims stated herein should be read as including those elements which are not necessary to the invention yet are in the prior art and are necessary to the overall function of that particular claim, and should be read as including, to the maximum extent permissible by law, known functional equivalents to the elements disclosed in the specification, even though those functional equivalents are not exhaustively detailed herein.

Although the present invention has been described chiefly in terms of the specific aspects of the invention, it is to be understood that the disclosure is not to be interpreted as limiting. Various alterations and modifications will no doubt become apparent to those skilled in the art after having read the above disclosure. Such modifications may involve other features which are already known in the design, manufacture and use of wireless electromagnetic communications networks, systems and MIMO networks and systems therefore, and which may be used instead of or in addition to features already described herein. The algorithms and equations herein are not limiting but instructive of aspects of the invention, and variations which are readily derived through programming or mathematical transformations which are standard or known to the appropriate art are not excluded by omission. Accordingly, it is intended that the appended claims are interpreted as covering all alterations and modifications as fall within the true spirit and scope of the invention in light of the prior art.

Additionally, although claims have been formulated in this application to particular combinations of elements, it should be understood that the scope of the disclosure of the present application also includes any single novel element or any novel combination of elements disclosed herein, either explicitly or implicitly, whether or not it relates to the same invention as presently claimed in any claim and whether or not it mitigates any or all of the same technical problems as does the present invention. The applicants hereby give notice that new claims may be formulated to such features and/or combinations of such features during the prosecution of the present application or of any further application derived therefrom. 

1. A method, comprising: channelizing a received navigation signal into a plurality of subband signals, each subband signal comprising a plurality of frequency channels; for each subband signal, computing linear combiner weights for the plurality of frequency channels based on one or more exploitable symbol stream properties of the received navigation signal; using the linear combiner weights to combine the frequency channels and excise interference in the each subband signal, thereby increasing signal-to-noise-and-interference of the each subband signal; and combining the plurality of subband signals to produce at least one of a detection statistic and a geo-observable estimate of the received navigation signal.
 2. The method of claim 1, wherein channelizing employs an analog-to-digital converter (ADC) configured to perform symbol-rate synchronous reception and channelization of the received navigation signal.
 3. The method of claim 1, wherein channelizing employs an analog-to-digital converter sampling rate that is equal to an integer multiple of the navigation signal’s baseband symbol rate.
 4. The method of claim 1, wherein channelizing employs an analog-to-digital converter sampling rate that is less than the received navigation signal’s bandwidth.
 5. The method of claim 1, wherein channelizing employs a fast Fourier transform (FFT), and FFT output bins are employed as the plurality of frequency channels.
 6. The method of claim 1, wherein an interference-excising linear algebraic combiner is configured to perform at least one of using the linear combiner weights to combine the frequency channels and combining the plurality of subband signals.
 7. The method of claim 1, wherein channelizing comprises varying the plurality of frequency channels in order to vary a number of processing degrees of freedom in the each subband signal.
 8. The method of claim 1, wherein the plurality of frequency channels is selected to be equal to or greater than a number of interferers in the each subband signal.
 9. The method of claim 1, wherein the one or more exploitable symbol stream properties comprises at least one of a periodic signal, known content in the symbol stream, a known pilot, a known modulation property, and a constant-modulus structure.
 10. The method of claim 1, wherein the plurality of subband signals are contiguous or separated in frequency, and the plurality of frequency channels in the each subband signal are contiguous or separated in frequency.
 11. The method of claim 1, further comprising using known positions of transmitters that transmit the at least one navigation signal in order to determine at least one of clock rate offset and ephemeris of a receiver that performs the method.
 12. The method of claim 1, further comprising estimating geo-observables from the interference, and geolocating one or more sources of the interference from the geo-observables.
 13. The method of claim 1, further comprising communicatively coupling to a data service or at least one navigation signal receiver for receiving third-party baseband symbol data or navigation data.
 14. The method of claim 1, further comprising computing positioning, navigation, and timing (PNT) analytics using the at least one geo-observable estimate, wherein the PNT analytics comprises at least one of blind Resilient PNT (RPNT) analytics, non-blind RPNT analytics, pilot-exploiting RPNT, rate synchronization, geolocation, navigation signal geo-observable estimates, navigation signal quality estimates, accuracy of the geo-observable estimates, and accuracy of the navigation signal quality estimates.
 15. The method of claim 14, wherein the computing performs analytics over integration times that are longer than a single baseband symbol, navigation symbol, or pilot period. 